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This  report  covers  Phase  I  of  a  study  of  Large  Waves  in  Channels, 
sponsored  by  DNA  under  the  SBIR  Program.  The  first  phase  was 
intended  to  review  the  knowledge  of  the  behavior  of  very  large 
waves  in  natural  channelized  geometries.  The  waves  considered 
are  very  long,  as  might  be  generated  by  submerged  explosions  of 
moderate  yield.  Owing  to  channelization,  such  waves  may  propa¬ 
gate  considerable  distances  with  little  attenuation.  On  sloping 
banks,  the  waves  are  subject  to  a  substantial  transverse  amplifi¬ 
cation,  further  increasing  their  potential  effects.  Breaking  may 
occur  in  the  bank  regions,  subjecting  shoreline  structures  and 
facilities  to  large,  destructive  effects. 


Having  demonstrated  the  potential  hazard  of  channel  waves,  it  is 
recommended  that  Phase  II  further  pursue  the  three-dimensional 
wave  theory  in  order  to  better  define  waye  transformation  ,  to 
remove  certain  theoretical  limitations,  to  develop  expressions 
for  water  particle  velocities  and  accelerations  in  the  bank 
regions  (for  assessment  of  forces  on  structures),  and  to  inves¬ 
tigate  the  conditions  under  which  the  bank  waves  will  break, 
forming  turbulent  bores. 
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1.  INTRODUCTION 


In  this  Phase  I  report  we  present  a  preliminary  review  of  the 
behavior  of  'Large  Waves  in  Channels.'  The  waves  of  interest  are 
long  oscillatory  cr  solitary  waves  such  as  might  be  generated  by 
large  explosions,  landslides,  or  other  impulsive  mechanisms.  We 
do  not  consider  waves  of  extreme  length  such  as  tides  or  flood 
waves,  nor  short  waves  in  the  range  of  the  normal  oceanographic 
environment.  The  waves  of  interest  are  those  that  are  longer 
than  the  width  of  the  channel  and  long  with  respect  to  the  depth, 
but  not  'too'  long. 

The  purpose  of  this  small  Phase  I  study  is  to  assess  the  general 
nature  of  this  problem,  to  review  available  information,  and  to 
define  directions  for  new  research  to  be  undertaken  in  Phase  II; 
Phase  I,  itself,  is  not  considered  a  research  phase.  Neverthe¬ 
less,  we  have  felt  it  appropriate  to  present  some  new  results 
even  here,  as  will  be  shown  later. 


Lituya  Bay  is  situated  on  the  Gulf  of  Alaska  south  of  Yakutat. 
As  shown  in  Figure  1,  it  is  roughly  T-shaped,  approximately  eight 
miles  long,  and  two  miles  wide;  depths  of  several  hundred  feet 
are  common.  The  landslide  (documented  in  Miller,  1960)  took 
place  in  the  northeast  portion  of  Gilbert  Inlet  (shown  shaded  in 
Figure  1).  Approximately  40  million  cubic  yards  of  rock  des¬ 
cended  at  high  speed  into  the  Bay  with  devastating  results. 

The  southwest  bank  of  Gilbert  Inlet  directly  opposite  the  slide 
was  inundated  by  the  slide-generated  water  waves,  to  a  height  of 
over  1500  feet,  completely  leveling  forest  to  that  elevation 
(Iida,  et  al,  1967).  This  is  certsinly  dramatic,  but  perhaps  not 
unexpected  considering  the  proximity  of  the  landslide  —  the 
opposite  bank  was  in  the  direct  path  of  the  'splash'  generated  by 
the  intruding  mass. 

More  interesting  were  the  effects  outside  the  immediate  region  of 
the  slide,  along  the  banks  of  the  central  portion  of  the  Bay. 
According  to  Iida,  et  al  (1967)  and  Wiegel  (1964),  a  large  wave 
progressed  seaward  through  the  Bay  with  heights  of  200-300  feet, 
as  estimated  from  the  elevation  to  which  large  trees  "were 
removed  as  if  they  were  matchsticks"  from  the  side  slopes 
(Wiegel,  1964). 

The  dotted  contours  in  Figure  la  show  the  elevation  of  destruc¬ 
tion,  the  numeric  labels  being  height  in  feet.  On  the  northern 
bank  one  sees  heights  of  about  90  to  200  feet,  while  on  the 
southern  bank  damage  elevations  range  from  about  100  to  600  feet. 
The  Bay  both  before  and  after  the  slide  is  shown  in  Figure  lb. 
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reproduced  from  Wiegel  (1964).  These  photos  show  the  large  areas 
of  forest  destroyed  by  the  passing  wave. 


There  are  two  major  reasons  why  waves  in  channels  pose  a  special 
problem.  First,  the  channelized  geometry  conducts  the  wave 
energy  from  the  source  to  a  site  a  distance  away  with  very  little 
loss.  This  is  in  contrast  to  the  usual  case  of  explosion  waves 
in  open  water.  In  open  water,  an  impulsive  wave  system  is  free 
to  expand  in  all  directions;  radial  spreading  causes  a  wave 
height  reduction  which  is  proportional  to  the  square-root  of 
distance,  leading  to  rapid  decay  of  height.  As  ordinary  waves 
approach  a  slope  they  refract  and  their  heights  are  limited  by 
breaking.  But  in  a  channel,  the  causes  of  height  reduction  are  I 

limited  to  dispersive  spreading  (small  for  long  waves)  and  energy 
loss.  The  major  source  of  energy  loss  may  be  turbulent  breaking  ’ 

in  the  zones  of  shallow  flow  on  the  banks.  But  these  breaking 

zones  are  precisely  the  damage  zones,  so  that  breaking  should 
actually  be  considered  a  damage  mechanism,  not  a  mitigating 
damping  mechanism.  The  wave  height  in  the  bank  zone  is  replen¬ 
ished  by  lateral  flow.  These  differences  are  summarized  in  the 
schematic  sketches  of  Figure  3. 

The  second  factor  contributing  to  a  high  hazard  in  channels  is 
the  possibility  of  significant  amplification  of  wave  height  on 
the  banks.  As  we  shall  show  later,  our  Phase  I  work  indicates 
that  amplifications  of  a  factor  of  two  or  more  are  to  be  expected 
on  typical  bank  slopes  for  moderately  wide  channels.  This 
amplification  is  found  as  a  higher-order  approximation  to 
solitary  wave  theory.  The  three-dimensional  solitary  wave  theory 
described  later  reveals  that  a  large  transverse  variation  of  wave 
amplitude  is  characteristic  of  solitary  waves  in  channels  and 
that  the  degree  of  amplification  at  the  banks  depends  upon  both 
channel  width  and  bank  slope.  The  general  nature  of  the  trans¬ 
verse  variation  is  found  to  be  consistent  with  the  limited 
available  experimental  data. 
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Figure  3b:  Wave  amplitude  variation  for  the  cases  shown  in 
Figure  3a. 
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The  explosion  (or  landslide)  waves  discussed  here  are,  in  a 
sense,  a  more  focused  phenomenon.  Under  proper  circumstances, 
they  will  propagate  with  almost  no  reduction  in  height  for 
considerable  distances.  Their  effects  are  greatest  along  the 
shorelines  (the  banks),  rather  than  in  deeper  water.  The  major 
dissipation  is  the  destructive  near-shore  breaking  which,  being 
of  limited  extent  (concentrated  right  at  the  bank),  may  not 
attenuate  the  wave  'too'  fast. 


The  major  content  of  this  report  is  contained  in  Sections  3  and  4 
where  we  describe  and  implement  a  three-dimensional  solitary  wave 
theory  to  estimate  the  lateral  enhancement  of  wave  height.  We 
have  worked  out  the  analytical  solution  for  triangular  channel 
cross-sections  and  have  constructed  a  simple  numerical  model  for 
the  simulation  of  a  range  of  trapezoidal  channel  geometries.  It 
is  found  that  amplification  factors  of  two  or  more  are  realistic 
for  natural  channels.  We  also  briefly  describe  the  effects  of 
friction  and  channel  bifurcation  as  dissipative  mechanisms. 
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Our  Phase  I  study  shows  the  importance  of  channel  effects  and 
indicates  the  direction  of  work  required  for  a  better 
understanding  of  this  sort  of  wave  behavior.  In  particular,  no 
investigation  of  particle  velocities  in  the  bank  zones  has  been 
undertaken,  and  no  information  is  available  describing  the 
breaking  characteristics  bn  the  banks.  Both  of  these  *questions 
are  vital  for  an  understanding  of  dynamic  effects.  There  is  also 
almost  no  experimental  data  to  guide  the  development  of  theories 
or  to  verify  them.  Consequently,  we  conclude  by  recommending  a 
program  of  theoretical  and  numerical  work  supplemented  by  small 
scale  laboratory  tests  for  a  Phase  II  continuation. 
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2.  REVIEW  OF  PRIOR  WORK 


2.1.  Source  Mechanisms 

We  consider  two  primary  source  mechanisms  for  the  sort  of  waves 
under  consideration  here:  explosions  and  landslides.  In  both, 
active  research  is  underway  elsewhere  and  has  reached  asstate  of 
adequate  understanding  for  the  purposes  of  this  study.  Conse¬ 
quently,  our  concern  is  limited  to  adopting  appropriate  order-of- 
magnitude  data  as  input. 

2.1.1.  Explosion  Sources 

Explosion  wave  prediction  (see  LeMehaute,  1971,  for  a  general 
summary)  is  accomplished  by  calibration  of  idealized  models  of 
wave  generation  by  an  initial  free-surface  deformation.  The 
source  is  usually  assumed  to  be  a  symmetric  'crater*  which 
collapses  under  gravity  yielding  an  outgoing  wave  system. 

The  initial  crater  is  fictitious  in  that  its  shape  is  chosen  so 
that  the  mathematics  is  tractable,  and  its  dimensions  are  chosen 
so  that  predicted  wave  heights  and  per iods  “match  observations 
reasonably  well.  The  predicted  waves  are  obtained  by  the 
transform  methods  originally  introduced  by  Kranzer  and  Keller 
(1959). 

*  v  v 

Considerable  effort  has  led  to  the  establishment  of  reasonable 
shapes  for  the  source  crater  as  well  as  to  empirical  scaling  laws 
relating  crater  dimensions  to  weapon  yield  and  placement.  For 
example,  in  deep  water,  the  familiar  curve  shown  in  Figure  4  is 
obtained,  relating  the  product  of  maximum  wave  amplitude,  n  ,  and 
range,  r  (this  product  being  conserved  in  deep  water)  to  ex¬ 
plosion  yield,  W,  and  depth  of  burst,  z.  Although  this  data 
shows  a  narrow  peak  in  wave  height  for  shallow  bursts  (the  so- 
called  upper  critical  depth  effect)  it  is  adequate  for  our 
purposes  to  follow  LeMehaute  (1985)  and  adopt  an  average  value 
given  by 


Hm ax  -  8W0.54/r 

in  which  the  dimensions  are  feet  and  pounds  of  TNT-equi valent . 
Similarly,  one  can  take  for  the  wave  period  T  (in  seconds) 

T  -  1.7W0.15 

from  which  the  length  can  be  calculated. 

These  deepwater  relationships  will  not  generally  be  appropriate 
for  waves  in  channels,  however,  owing  to  the  relative  shallowness 
of  most  real  channels.  The  effect  of  shallow  water  is  to  greatly 
reduce  the  efficiency  of  wave  generation. 

Wang,  et  al  (1977)  reviewed  available  explosion  data  for1  shallow 
water  (especially  the  data  of  Strange,  1955)  and  revised  the 
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scaling  laws  accordingly.  As  given  in  the  recent  summary  of 
LeMehaute  (1985),  Wang's  revised  expressions  are 

8  0.25  +  8/3  1/3  0.93 

Hmax1- /( W  )  -  1.44  (d/  W  ) 

1/3  0.007 

8  =  0.83  (d/  W  ) 

f  1/3 

These  expressions  are  appropriate  for  reduced  depths  ( d /  W  ) 
less  than  one  (the  original  data  being  from  0.088  to  0.585). 

The  deepwater  relations  given  above  can  be  accepted  for  reduced 
depths  greater  than  about  six.  For  intermediate  depths, 
LeMehaute  (1971)  suggests  that  heights  should  be  reduced  by  half 
between  a  reduced  depth  of  6.0  and  a  depth  of  1.0;  a  simple 
linear  variation  is  suggested. 

An  important  question  for  channel  waves  is  the  dependence  of  wave 
length  on  water  depth  for  a  given  explosion  yield.  In  order  to 
exist  and  survive,  these  waves  must  be  longer  than  both  the 
channel  depth  and  the  channel  width. 

It  seems  reasonable  on  physical  grounds  to  argue  that  the  wave 
length  of  the  highest  waves  generated  by  an  explosion  will  be,  if 
not  independent  of,  at  least  insensitive  to  the  water  depth.  The 
dominant  length  scale  in  the  spectrum  of  waves  generated  by  an 
initial  'crater'  must  be  of  the  same  order  of  magnitude  as  a 
characteristic  horizontal  measure  of  the  crater  itself.  That  is, 
we  should  expect  energy  to  be  concentrated  at  lengths  near  the 
diameter  of  the  crater,  within  a  small  factor. 

Indeed,  data  shows  that  in  deep  water  the  length  of  the  largest 
waves  is  about  three-fourths  the  assumed  crater  diameter, 
confirming  our  assumption. 

Now,  in  shallow  water  it  is  reasonable  to  assume  that  to  a  first 
approximation  the  horizontal  scale  of  the  crater  is  unaffected  by 
the  depth.  That  is,  the  blast  creates  a  disturbance  which  has 
about  the  same  diameter  no  matter  what  the  depth  is,  although  the 
vertical  crater  scale  might  be  substantially  affected. 

In  the  absence  of  a  better  simple  alternative,  we  shall  then 
assume  that  wavelengths  in  channels  can  be  estimated  from 
deepwater  empirical  relationships.  Should  it  be  necessary  to 
estimate  wave  periods,  one  can  simply  work  backward  from  the 
length  and  the  depth. 

2.1.2  Landslides 

A  detailed  numerical  model  of  landslide  generated  water  waves  was 
developed  by  Chiang,  Divoky,  Parnicky,  and  Wier  (1981).  This 
model  incorporates  a  sophisticated  'box-car'  model  of  landslide 
dynamics  on  an  arbitrary  slope,  accounting  for  the  interaction  of 
the  slide  with  the  water  upon  entry  into  the  basin.  This  coupled 
calculation  then  proceeds  to  follow  the  generation  and  propaga- 
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tion  of  the  water  waves.  Other  models  of  varying  complexity  have 
been  developed  by  a  number  of  authors  including  especially  Raney 
and  Butler  (1975)  and  Noda  (1969),  the  latter  being  an  analytical 
development  describing  the  intrusion  of  a  mass  of  simple  shape 
into  a  fluid. 

It  had  been  hoped  that  the  model  of  Chiang,  et  al,fcould  be 
converted  from  mainframe  to  PC  for  the  present  study  and  that 
simple  example  cases  could  be  run  for  an  idealized  geometry. 
However,  this  proved  impossible  within  the  limited  Phase  I 
resources  . 

More  than  adequate  for  our  purposes,  however,  are  empirical 
guidelines.  Slingerland  and  Voight  (1982)  present  a  simple 
regression  model  relating  landslide  induced  wave  height  to  slide 
kinetic  energy  and  water  depth.  Their  expression  is 

l°g(nmax/d)  "  a  +  b  1 °8 C E ) 

in  which  a  and  b  are  regression  coefficients  and  E  is  a  dimen¬ 
sionless  slide  kinetic  energy  given  hy 

E  -  1/2  MV2/(Pgd4) 

Here,  M  is  the  slide  mass,  g  is  the  acceleration  of  gravity,  and 
V  is  the  maximum  slide  velocity.  The  wave 'amplitude  rimax  is  the 
value  observed  at  a  standard  reference  distance  of  4d  from  the 
slide  point  of  entry,  and  assumes  that  the  waves  spreads  freely 
into  a  180-degree  half-space.  Should  the  slide  occur  at  a 
corner,  Slingerland  and  Voight  recommend  the  artifice  of  simply 
adjusting  the  assumed  energy  in  order  to  account  for  the  effects 
of  reduced  spreading  area. 

Figure  5  shows  the  regression  relation  along  with  the  data  points 
upon  which  it  is  based.  It  is  clear  that  within  this  data  base, 
the  simple  relation  suggested  by  Slingerland  and  Voight  is 
entirely  adequate. 

To  use  this  formula  for  a  first  order  estimate  in  channels  (with 
slide  intrusion  assumed  to  occur  on  one  side),  the  amplitude 
should  be  adjusted  to  account  for  channel  width.  It  is  reason¬ 
able  to  imagine  that  the  height  which  should  be  considered  the 
initial  height  for  propagation  in  a  channel  is  approximately  that 
obtained  when  the  expanding  semi-circular  wave  encounters  the 
opposite  bank.  At  this  time,  complicated  reflection  and  dissipa¬ 
tion  processes  on  the  opposite  bank  truncate  the  wave  expansion 
and  leave  two  segments  of  the  wave  to  propagate  up-  and  down¬ 
stream.  These  segments  are  perpendicular  to  the  bank  on  the  near 
shore  but  are  curved  toward  the  bank  on  the  opposite  shore, 
leading  to  additional  dissipation.  However,  if  the  wave  is 
reasonably  long  compared  to  the  channel  width  it  may  be  assumed 
that  it  will  tend  toward  cross-channel  uniformity  with  distance. 

Since  the  regression  expression  of  Slingerland  and  Voight  adopts 
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a  reference  distance  of  4d ,  the  height  prediction  must  be  scaled 
by  a  factor  such  as 


(4d/B0)P 

in  which  B0  is  the  channel  width  and  p  is  an  exponent  ranging 
between  0.5  (longwaves,  no  dispersion)  and  1.0  (fully  (fispersive 
waves).  It  is  recommended  that  p  be  taken  as  0.5  unless  the 
channel  is  quite  deep;  this  choice  will  tend  to  slightly  underes¬ 
timate  the  wave  height,  helping  to  compensate  for  unconsidered 
dissipation . 
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2.2.  Wave  Behavior 


2.2.1.  Theoretical  Descriptions 

Very  little  of  the  literature  on  waves  in  channels  is  pertinent 
to  this  study,  most  dealing,  for  example,  with  the  propagation  of 
flood  waves  or  dam-break  waves.  These  are  phenomena  quite 
distinct  from  the  explosion-scale  intermediate  waves  of  interest 
her  e . 

However,  a  small  number  of  papers  were  found  which  are  directly 
applicable,  developing  a  theory  of  cnoidal  and  solitary  waves  in 
channels  of  arbitrary  cross-section.  When  taken  to  higher  order, 
these  theories  are  capable  of  describing  major  features  of  flow 
in  the  sloping  bank  regions  which  are  the  areas  of  primary 
interest  in  this  study  (up  to  the  point  of  breaking,  that  is). 

The  papers  of  interest  followed  the  important  earlier  work  of 
Korteweg  and  de  Vries  (1895),  Ursell  (1953),  and  Benjamin  and 
Lighthill  (1954)  who  laid  the  necessary  groundwork.  Korteweg  and 
de  Vries  and  Ursell  had  clarified  the  roles  played  by  the  two 
fundamental  dimensionless  length  scales  involved  in  longwave 
theory.  These  are  a  and  e  defined  by 

o  «  h/L  e  -  H/h 

•  '  v 

in  which  h  is  a  typical  depth,  L  is  a  horizontal  length  scale 
such  as  wavelength,  and  n  is  a  typical  amplitude. 

Ursell  showed  that  the  simultaneous  smallness  of  o  and  £  does  not 
lead  to  a  unique  theory  as  was  assumed  (the  longwave  paradox), 
but  that  the  relative  magnitudes  of  these  parameters  is  equally 
important.  The  combined  parameter 

U  -  e/O2 

is  called  the  Ursell  parameter  and  can  be  interpreted  as  a 
measure  of  the  relative  importance  of  nonlinearity,  which  tends 
to  steepen  the  wave  with  propagation  distance  (measured  by  e  ) , 
and  frequency  dispersion,  which  tends  to  broaden  the  wave 
(measured  by  o  ) .  It  was  shown  that  when  the  parameter  U  is  on 
the  order  of  unity,  these  opposing  tendencies  may  balance, 
permitting  the  existence  of  waves  of  permanent  form. 

Korteweg  and  de  Vries  (1895)  and  Benjamin  and  Lighthill  (1954) 
developed  the  governing  equations  for  cnoidal  waves  (and,  in  the 
limit,  solitary  waves)  ,  the  latter  in  terms  of  gross  flow 
quantities . 

Peters  (1966)  considered  solitary  waves  in  a  channel  of  arbitrary 
cross-section.  His  paper  is  of  interest  since  he  first  develops 
the  rotational  case  and  then  the  irrotational  case  as  a  limit. 
He  also  gives  an  analytical  solution  for  semi-circular (channels 
which  will  be  discussed  in  a  later  section. 


The  papers  of  most  interest  are  those  of  Peregrine  (1968)  and 
Fenton  (1973).  The  work  of  these  authors  is  quite  similar.  Both 
consider  long  waves  in  arbitrary  channels  and  derive  similar 
expressions  for  flow  quantities.  Peregrine's  development  is  the 
more  standard,  being  in  terms  of  particle  velocities  and  eleva¬ 
tion  whereas  Fenton's  approach  is  really  an  extension  of  that  of 
Benjamin  and  Lighthill,  being  defined  in  terms  of  integrated  flow 
quantities . 

Although  we  shall  show  Peregrine's  development  in  a  later  section 
(upon  which  we  have  based  a  preliminary  model  of  transverse  wave 
amplitude  variation),  Fenton's  derivations  present  certain 
conceptual  advantages.  In  particular,  he  retains  both  o  and  £  in 
his  derivation,  thereby  showing  precisely  their  relative  impor¬ 
tance.  Peregrine,  on  the  other  hand,  simply  sets  U  ■  1  and 
suppresses  o  in  favor  of  £  . 

These  papers  do  not  go  much  beyond  development  of  the  governing 
equations,  except  to  discuss  the  transverse  amplitude  variation 
in  a  general  way.  No  pursuit  of  particle  velocities  is  made,  for 
example,  such  as  will  be  necessary  to  define  breaking  criteria 
and  dynamic  effects. 

2.2.2.  Observations 

*  '  v 

Experimental  observations  of  the  important  features  of  waves  in 
channels  are  almost  totally  lacking.  In  particular,  we  are 
interested  in  the  transverse  wave  profile  and  the  likelihood  that 
amplitudes  on  the  banks  may  be  significantly  greater  than  at  the 
center  of  the  channel. 

It  is  not  clear  why  this  has  not  been  a  topic  of  interest.  A 
large  quantity  of  experimental  data  for  trapezoidal  channels  is 
reported  in  Benet  and  Cunge  (1971)  which  was  expected  to  be  a 
good  source  of  data  for  this  review.  Surprisingly,  however,  the 
critical  data  is  not  reported.  All  measurements  seem  to  have 
been  made  on  the  banks  and  what  appears  to  be  amplification  data 
is  actually  the  ratio  of  peak  bank  height  to  the  surcharge  in 
mean  water  level  associated  with  a  sudden  discharge. 

Data  for  waves  impinging  at  a  small  angle  to  a  slope  was  obtained 
by  Divoky  and  Lane  (1973),  but  is  of  limited  use  here.  Firstly, 
the  waves  were  not  long  with  respect  to  the  depth  and  secondly, 
the  test  arrangement  was  in  a  2-D  basin,  not  a  channel.  Conse¬ 
quently,  the  waves  in  many  cases  were  not  propagating  stably,  the 
observed  amplifications  showing  considerable  scatter. 

Nevertheless,  it  may  be  of  interest  to  show  the  most  pertinent 
cases.  Figures  6  and  7  show  the  observed  amplitude  profiles 
along  a  crest  for  nominal  angles  of  incidence  of  0  and  3.5 
degrees,  for  three  nominal  beach  slopes:  20,  30,  and  45  degrees. 
The  water  depth  at  the  toe  of  the  slope  was  1.2  ft  and  rthe  wave 
period  was  about  1.1  seconds.  Amplitudes  along  the  crests  were 
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measured  with  a  resistance  gage  or,  on  the  slope,  estimated  from 
the  limit  of  the  wetted-zone.  Amplitudes  were  arbitrarily 
normalized  to  unity  at  the  farthest  offshore  points. 

It  can  be  seen  from  these  figures  that  for  both  30  and  45  degrees 
there  is  an  amplitude  enhancement  of  as  much  as  a  factor  of  four. 
On  the  20  degree  slopes,  however,  amplifications  were  only  about 
half  as  large.  The  reasons  for  this  were  not  evident,  especially 
since  the  same  trend  is  observed  for  both  breaking  and  non¬ 
breaking  waves,  but  it  does  seem  to  be  a  consistent  feature  of 
the  data  for  waves  of  grazing  incidence.  Perhaps,  because  of  the 
small  experimental  scale,  scale  effects  became  disproportionate 
in  the  shallow  wedge  region  of  the  bank  in  the  20  degree  case. 

The  wave  generator  was  adjusted  to  produce  both  small  waves 
showing  no  near-shore  breaking,  as  well  as  larger  waves  revealing 
breaking  at  the  tip  of  flow.  It  was  found  that  for  0  degree 
incidence,  wave  breaking  occurred  for  a  deepwater  wave  steepness 
H0/L0  of  0.01,  0.013,  and  0.017  for  beach  slopes  of  20,  30,  and 
45  degrees,  respectively.  The  same  trend  was  observed  for  an 
incident  angle  of  3.5  degrees,  with  each  critical  steepness  about 
0.003  higher  than  for  the  0  degree  case.*  For  less  oblique 
incidence,  the  breaking  criterion  became  extremely  erratic 
reflecting  the  complex  behavior  in  a  Mach-stem  zone. 

The  final  --  and  most  pertinent  --  data'- is  that  obtained  by 
Peregrine  (1969)  in  two  small  lab  test  series.  The  results  are 
semi-quantitative,  being  presented  only  as  two  drawings  of 
surface  profile  across  half  trapezoids  for  various  center 
amplitudes.  Both  cases  involved  side  slopes  of  45  degrees,  one 
having  a  bottom  half-width  of  1.5,  the  other  2.0  (normalized  by 
depth).  These  are  reproduced  in  Figure  8.  The  elevations  shown 
at  the  left  edge  (full-channel  center  line)  are  normalized  by 
depth.  The  solid  lines  represent  the  crest  profiles,  while  the 
small  x's  indicate  the  corresponding  values  from  Peregrine's 
model.  In  general,  the  agreement  is  seen  to  be  very  good. 
Surface  heights  at  the  banks  appear  to  be  about  50%  greater  than 
at  the  channel  center  for  the  larger  amplitudes. 

Note  in  particular  that  in  the  case  of  a  bottom  width  of  2.0  and 
a  center  amplitude  of  0.27,  the  crest  at  the  bank  is  breaking,  or 
nearly  so.  Despite  this,  the  theoretical  prediction  appears 
quite  good. 
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Figure  8a:  Transverse  Profiles  from  Peregrine,  1969 
Bottom  Width  =  1.5,  Slope  »  45°. 
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Transverse  Profiles  From  Peregri 
Bottom  Width  ■  2.0,  Slope  ■  45 
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3.  WAVES  IN  A  UNIFORM  CHANNEL 

3.1  Nonlinear  Three-dimensional  Waves 


3.1.1  Theoretical  Development 

We  have  adopted  the  development  given  by  Peregrine  (1968)  for  use 
in  this  study.  As  noted  earlier.  Peregrine's  formulation  is  more 
immediately  useful  in  our  context  than  is  the  presentation  of 
Fenton,  for  example,  although  the  essential  features  remain  the 
same.  Consequently,  this  summary  follows  Peregrine  closely. 


Consider  a  cartesian  coordinate  system  embedded  in  the  still 
water  surface  with  the  x-axis  taken  in  the  direction  of  wave 
propagation,  the  y-axis  taken  transverse  to  the  channel,  and  the 
z-axis  positive  upward  (Figure  9  ).  The  y-origin  may  be  visual¬ 
ized  as  above  the  channel  invert  or,  in  our  later  consideration 
of  symmetrical  channels,  at  the  midpoint  of  the  cross-section. 
We  consider  the  flow  to  occur  in  a  channel  of  arbitrary,  but 
uniform,  cross-section;  that  is,  the. cross-section  is  defined  by 

zc  ■  zc(y)  ( i) 


which  is  independent  of  x. 

•  '  v 

Introduce  dimensionless  variables 

(x.y.z)  -  (x*,y*,z*)/h0 
t  -  t*  /(g/h0) 

(u.v.w)  -  u  -  u*  /  /(gh0) 

P  -  (P*  -  P0)/P8h0  (2) 


in  which  the  starred  quantities  are  dimensional,  g  is  the 
acceleration  of  gravity,  h0  is  a  characteristic  water  depth  (such 
as  the  maximum  undisturbed  depth),  (u,v,w)  are  the  components  of 
the  fluid  particle  velocity  H ,  p  is  pressure,  and  p0  is  at¬ 
mospheric  pressure. 


Neglecting  only  the  shear  stresses,  the  momentum  equations  become 
simply  the  Euler  equations, 

(ft  +  (fr.V)if+  Vp  +  (0,0,1)  -  0  (3) 


and  the  incompressible  continuity  equation  is  just 

V*  u  -  0  (4) 


This  can  also  be  written  in  terms  of  the  macroscopic  quantities  Q 
and  A  (volume  flux  and  instantaneous  wetted  cross-section)  as 

A t  +  Qx  "0  ( 4a ) 
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Figure  9:  Sketch  of  Channel  Coordinate  System. 
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We  obviously  have 


A  »  A(x,t)  and  Q  =  u  dydz 

We  assume  irrotational  flow  so  that 

V  Xu  =  0  (5) 

The  free  surface  is  denoted  by  n;  on  the  surface 

z  =  n(x,y,t) 

where  we  must  satisfy  the  boundary  conditions 

p  a  0 

Ht  +  UX  +  Vy  a  W  (6) 

Note  that  in  these  equations  and  those  to  follow,  we  adopt  the 
subscript  convention  to  denote  partial  differentiation  (in  the 
interests  of  notational  simplicity). 

Two  dimensionless  relative  length  scales  are  now  introduced  in 
order  to  proceed  with  an  expansion  of  the$£  equations.  These  are 
the  relative  depth,  o  ,  and  the  relative  amplitude, e  ,  given  by 

0-  h0/L* 

e -  a*/hQ  (7) 

in  which  L*  is  a  characteristic  wavelength  and  a*  is  a  measure  of 
the  peak  amplitude. 

In  our  problem,  we  assume  that  o  is  very  small  since  we  are 
restricting  attention  to  long  waves.  We  will  also  assume  that  e 
is  small,  although  we  shall  be  interested  in  waves  which  may  be  a 
few  tenths  of  the  depth  in  height. 

A  number  of  authors  have  shown  (for  constant  depth  or  rectangular 
channels)  that  solitary  waves  of  unchanging  form  are  possible 
only  if  a  balance  exists  between  two  processes  measured  by  o  and 
c.  The  effect  of  finite  amplitude,  measured  by  e,  is  to  cause  a 
steepening  of  wave  form  with  propagation  distance,  eventually 
leading  to  breaking.  The  effect  of  frequency  dispersion, 
measured  by  o,  is  the  opposite,  tending  to  reduce  average 
steepness  with  distance.  If  these  opposing  influences  are 
comparable,  then  a  finite-amplitude  wave  of  permanent  form  is 
possible.  Not  only  solitary,  but  also  periodic  waves  --  the 
cnoidal  wave  approximation  --  are  found  in  this  case,  although  we 
shall  in  this  preliminary  review  focus  on  the  solitary  case. 

Ursell  (1953),  in  particular,  clarified  the  importance  of  the 
relative  magnitudes  of  a  and  e,  showing  that  the  ratio 
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e/  o  2 


a*L*2  /  h03 


(8) 


is  of  fundamental  importance.  A  unique  theory  is  not  guaranteed 
by  the  simultaneous  smallness  of  a  and  e.  Instead,  the  magnitude 
of  U  (the  Ursell  parameter)  must  also  be  specified.  Only  when  U 
is  of  order  unity  do  finite-amplitude  permanent  forms  arise. 

In  order,  then,  to  obtain  the  wanted  solution  Peregrine  adopts 
the  assumption  that  U  ■  1,  replacing  o  by  e  in  his  development. 
It  might  be  noted  that  in  the  similar  work  of  Fenton  (1973)  U  is 
also  assumed  to  be  of  order  unity  but  both  o  and  e  are  retained 
in  the  expansions,  greatly  clarifying  their  roles  and  the  order 
of  approximation  at  each  stage.  While  we  follow  Peregrine  here 
for  simplicity,  future  work  would  more  likely  be  based  on  the 
approach  of  Fenton. 

The  independent  variables  are  first  scaled  as 

xi  ■  x/e  tj  ■  t/e  (9) 

which  are  simply  x*/L*  and  t*/T,  where  T  is  a  characteristic  time 
scale  such  as  the  characteristic  length  L*  divided  by  the  speed 
of  longwave  propagation.  The  variables  n  and  u  vanish  in  the 
absence  of  waves  (neglecting  a  mean  current')  and  so  are  expanded 
in  the  form 


cf  \  +  e*f  2  +  •  •  . 

The  cross-sectional  velocities  are  expanded  as 

/e  (  ef  i  +  e*f  2  +  •  •  • ) 

and  the  variables  p  and  A,  not  vanishing,  are  taken  to  be  of  the 
form 


f0  ♦  Ef i  +  e*f 2  +  •  •  • 

The  expansion  solution  is  now  straightforward.  The  expanded 
forms  are  substituted  into  the  governing  equations,  terms  of  the 
same  order  are  grouped  and  equated  to  zero,  giving  successive 
approximations  to  the  dependent  variables. 

It  is  easy  to  see  that  the  only  zero’th  order  solutions  are 

p0  ■  -z  and  A0  ■  constant 

That  is,  the  pressure  is  hydrostatic  and  the  cross-sectional  area 
is  constant  (assuming  a  uniform  channel)  to  this  approximation. 

Peregrine  obtains  the  first  order  approximation  by  defining  a 
two-dimensional  velocity  potential  function  $1 ( *  1 , y , z ,x i )  such 
that 
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V1  -  *ly 


and 


(10) 


W1  =  *lz 

which  is  permitted  by  the  i r r o t a t i ona 1 i t y  condition  to  the  first 
order  ( Eq  5) 


v 1 z  -  w  1  y  =  0  (11) 

Equation  5  also  shows  that  to  first  order  u  is  independent  of  y 
and  z  so  that 


U1  =  ui(xi.ti) 

Using  Peregrine's  two-dimensional  gradient  operator  defined  as 

7i  ■  (  0  ,  3/3  y  ,  3/3  z  ) 

the  first  order  y  and  z  momentum  equations  yield 

7iPi  *  0 

Expanding  the  pressure  condition  on.  the  frq.e  surface  gives  the 
first  order  condition 


Po  +  EP1  ■  0  on  Z“£ni 

so  that  the  first  order  expression  for  p  is'  found  to  be 

Pl(xi.ti)  -  ni(xi.ti) 

The  x  momentum  equation  at  first  order  is  found  to  be 

(«i)tl  +  <ni)Xl  -  0  (12) 

and  the  integrated  form  of  the  continuity  equation  is,  similarly, 

(Al)tx  +  (Ql)xx  *  0 

Letting  the  width  of  the  channel  be  denoted  by  B(z),  it  is  clear 
that  the  first  correction  to  channel  cross-sectional  area  is 

Aj  -  B(0)nj  -  B0Hi 

Furthermore,  we  must  have  Qi  »  Ao  u  j  so  that  the  continuity 
equation  can  be  rewritten  as 

Bo(ni)tl  +  Ao  ( u  1 )  Xj  ■  0  (I3) 

Now,  uj  is  eliminated  from  Equations  12  and  13  giving 

(’llHjtj  *  A0/B0  ( H 1 ) Xj  (1^) 

This  is  just  the  wave  equation  from  which  we  can  identify  the 
first  approximation  to  the  wave  celerity,  c0,  to  be  given  by 
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(15) 


c0  *  A0/B0 

which  is  just  the  mean  depth  of  the  channel. 

Now  we  introduce  the  cross-sectional  velocity  potential  into  the 
continuity  equation  (Eq  4),  finding 

(ul)Xl  +  71  *1  -  o 

On  all  solid  boundaries  we  must  satisfy  (  <t>i)n  *  0  and  on  the 
surface  the  first  order  condition 

(*l>z  *  ~co  <ul)Xl  at  z  -  0 

where  an  H-derivative  has  been  replaced  by  a  u-derivative  using 
Equation  13. 

Guided  by  inspection  of  this  last  expression.  Peregrine  intro¬ 
duces  a  new  function,  <1>  ,  defined  by 

<t>(xi,y,z,ti)  -  -  (ui)x1'Ky.z) 

Comparison  with  the  foregoing  equations  shows  that  'l'  must  satisfy 

Vj  4<  ■  1  everywhere' 

vpn  —  0  on  solid  boundaries 

-  Cq  on  z  -  0  (16) 

Clearly,  this  is  a  well-posed  Neumann  problem  which  we  solve 
numerically  in  a  subsequent  section. 

To  obtain  the  transverse  variations,  it  is  necessary  to  proceed 
to  the  next  approximation.  Collecting  terms  through  the  second 
order  in  e  in  the  momentum  equation  gives 

<Vi*i)t  +  V1P2-O  (17) 

which  can  be  integrated  to  yield 

P2  -  Cui)Xl  +  DUl.ti) 

in  which  D  is  an  arbitrary  function.  The  free  surface  pressure 
condition  at  this  order  is  P2  ■  n  2  so  that  the  expression  above 
immediately  gives  the  transverse  wave  amplitude  variation  as 

n2  “  (ui)xj  tj  *  (7.0)  +  D  (18) 

In  other  words,  the  variation  of  the  free  surface  is  like  'l'(y.O) 
which  is  obtained  as  discussed  above.  r 

From  the  irrotationality  condition,  taking  the  y  and  z  components 
to  0(e2  ),  Peregrine  shows  the  next  approximation  for  u  to  be 
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u 2  *  -(ul)xl  Xj  ^(y.z)  +  E (  x  i  ,  1 1  ) 

in  which  E  is  arbitrary  (obtained  by  an  integration  as  in 
Equation  18). 

Proceeding,  the  next  order  in  the  x-momentum  equation  (toe2  )  is 
found  to  be 

(  u  2  )  t  +  u  1  (  u  i )  x  +  (  p  2  )  x  =  0 
i  i  i 

Neglecting  the  cross-sectional  variations  in  both  u2  and  P2  this 
becomes  an  equation  relating  D  and  E 

Et  +  ui(ui)x  +  Dx^  =  0  (20) 

Now,  assuming  that  over  the  small  range  from  the  still  water 
surface  to  the  instantaneous  surface,  the  top  width  is  given  by 

B(z)  «  B0  +  zB^ 

and  incorporating  the  transverse  variation  of  Equation  18,  the 
second  order  term  in  the  expansion  of  flow  area  is 

A 2  -  Bini/2  +  B0D  +  Bo^bC*!)^  (21) 

Here,  Bi  denotes  the  slope  of  the  bank  at  z  •  0,  and  'J'g  is  the 
mean  value  of  <i>  on  the  surface  given  by 

H>B  *  1/Bo  Jb  *(y.0)dy  (22) 

Similarly,  the  second  order  term  in  the  flux,  Q,  is  shown  to  be 

Q2  -  B0u i  +  A0E  -  A0^a(ui)Xi  Xi  (23) 

in  which  ^A  is  the  mean  value  of  over  the  entire  cross-section 
given  by 

♦A  -  1  / A o  jjA  ^(y.z)dydz  (24) 

Peregrine  now  introduces  another  transformation  of  the  dependent 
variables,  introducing  quantities  which  are  correct  to  the  second 
order,  but  neglecting  the  cross-sectional  variations.  That  is, 
for  the  free  surface,  define 

2 

5  -  eni  +  cD 

while  for  the  longitudinal  flow  define 

u '  ■  eu  i  +  e*E 

With  these  definitions,  one  now  combines  £  times  Equation  12  with£ 
times  Equation  20,  reverts  to  un-scaled  x  and  t  coordinates,  and 
finds  for  the  momentum  equation 
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if 


ut 

+  u’ux 

+  ;x  -  o 

(25) 

The  continuity  equation 

i  s 

not  as 

simple,  being 

(  C  +  k;2 /2 ) t  +  ( CqU  ' 

+ 

;u’)x  + 

('i'B  ~  ^a)(u'  )xtt  *  0 

*(26) 

in  which  k  is  Bj/Bq. 


A  solution  of  these  equations  which  is  extremely  similar  to  the 
familiar  solitary  wave  is 

i  1 

v2 


sech^ 


x  -  ct 


a  (1  -  kc73) 
o 


C*  (K  -  ) 

o  B  A 


in  which  the  wave  celerity  is  given  by 

c  -  c0[l  +  (1  -  kcj/3)a/2c£] 


(27) 


(28) 


Peregrine  (1968)  concludes  at  this  point  without  presenting  an 
explicit  expression  for  the  transverse  amplitude  variation.  The 
missing  steps  are  easy  to  supply,  however.  In  Equation  18  we 
replace  the  time-derivative  of  velocity  by,  the  x-derivative  of  n 
using  Equation  12.  This  second  derivative  of  n  is  then  evaluated 
using  Equation  27.  Since  we  are  interested  only  in  the  crests, 
we  then  set  the  phase  position  to  zero  and  find  the  transverse 
variation  of  amplitude,  a,  to  be 


.  .  ♦  °o  Bo  (1  -  (29, 

2A0  (4,B  ■  V 

This  result  is  given  both  by  Fenton  (1973)  and  in  a  second  paper 
by  Peregrine  (  1969),  the  latter  without  derivation.  In  this 
expression,  a0  is  the  wave  amplitude  at  the  minimum  of  'J'(y.O) 
which  we  shall  later  construe  to  be  at  the  center  of  a  simple, 
symmetric  cross-section  (in  fact,  at  y  ■  0). 

We  should  note  before  concluding  this  section  that  two  limita¬ 
tions  are  inherent  in  the  theoretical  development  (beyond  the 
requirement  of  small  0  and  e  ).  First,  it  was  tacitly  assumed 
that  the  bank  slope  is  not  too  small.  This  arises  in  Equa¬ 
tion  21:  if  B}  is  large  owing  to  an  extremely  gentle  slope,  its 
contribution  is  no  longer  second  order,  violating  the  assumption. 

The  second  limitation  is  that  the  channel  should  not  be  too  wide 
compared  with  the  depth.  Peregrine  (1969)  shows  that  contrary  to 
expectation,  increasing  the  channel  width  eventually  causes  such 
a  large  transverse  amplitude  variation  that  the  assumption  of 
small  amplitude  must  be  violated.  r 

This  completes  the  theoretical  development  of  the  topic  to  the 
extent  needed  for  this  preliminary  Phase  I  effort.  We  shall 
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suggest  that  in  Phase  II  we  pursue  this  theory  along  three  major 
lines  : 

*)  Investigation  of  particle  velocities  in  the  bank  area  for 
the  determination  of  forces  on  bank  structures  t 

*)  Investigation  of  both  velocities  and  accelerations  at  the 
crest  in  the  bank  area  in  order  to  define  wave -breaking 
criteria  (both  kinematic  and  dynamic) 

*)  Generalization  to  remove  or  mitigate  the  errors  as¬ 
sociated  with  increasing  channel  width  and  decreasing  bank 
slope 


s 
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3.1.2  Analytical  Solutions 


Analytical  solutions  of  the  foregoing  equations  have  been 
obtained  only  in  a  small  number  of  simple  geometries.  For 
example,  in  a  rectangular  channel  (physically  equivalent  to  the 
unbounded  constant  depth  case)  we  have 

*  (y  ,z)  -  <z) 

=  Cq  *  1  at  z  >  0 

^z*0  at  z  *  -1 

with  which  (guessing  a  polynomial  solution)  we  immediately  find 

z  +  z  2/  2 

Isolevels  of  ^  are  horizontal  and  consequently  all  cross-section¬ 
al  flow  is  vertical  (perpendicular  to  the  isolevels).  Equa¬ 
tion  29  shows  that  the  surface  is  level  since.  ^  is  not  a  function 
of  y  (and  we  have  tacitly  taken  the  arbitrary  constant  in  the 
solution  for  to  be  such  that  <i>  is  zero  on  the  free  surface). 
This  is,  of  course,  a  trivial  case. 

Somewhat  more  interesting  is  the  result  found  by  Peters  (1966) 
for  a  semi-circular  cross-section.  The  expression  for  ^  is 
moderately  complex,  and  gives  the  following  approximation  for 
bank  amplification  of  wave  amplitude 

ab/a0  » 1  +  a0/6r 

In  this  expression,  ab  is  the  amplitude  at  the  shore  and  r  is  the 
circular  channel  radius.  It  is  noted  that  the  amplification  is 
quite  small.  Even  for  near  limit  wave  amplitude,  the  ratio  is 
only  about  1.1,  or  a  10Z  increase  near  the  banks.  Physically, 
this  seems  to  be  due  to  the  fact  that  the  banks  are  vertical  at 
the  still  water  level,  and  nearly  vertical  over  a  substantial 
depth. 

A  third  analytical  solution  of  more  practical  interest  can  be 
found  for  channels  of  triangular  cross-section.  Guessing  a 
polynomial  expression  and  noting  that  c£  is  1/2  for  a  triangle, 
it  is  easy  to  find 

^  -  z/2  +  (y2  +  z^)/4 

Again,  the  arbitrary  constant  is  chosen  so  that  'l'  is  zero  at  the 
origin.  We  have  A0  ■  B0/2  in  this  case  and  Bi  ■  B0.  The  sloping 
bank  is  defined  by 

y  -  (1  +  z)B0  f 

which  becomes  the  upper  limit  on  the  inner  integral  for  'i'A 
(Equation  24).  We  find,  after  some  algebra 
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*B  ■  B20  /  12 
*A  =  (Bo  -  3)/24 


so  that  Equation  29  reduces  to 


a 


ao  + 


5  y  2 /  (3  +  Bq) 


Of  greatest  interest  is  the  amplitude  at  the  bank  (y  =  B0)  which, 
when  normalized  by  the  central  amplitude,  is  the  amplification 
factor  given  by 


R  -  afc/a0 


1  +  a0 


5B02/(3  +  Bq) 


Note  that  B0  is  just  tan(a)  where  a  is  the  angle  between  the  bank 
and  the  vertical.  We  shall  denote  the  angle  between  the  bank  and 
the  horizontal  by  0  . 


The  following  Table  shows  the  variation  of  R-with  both  side  slope 
and  center-channel  amplitude  for  triangular  channels. 

R 


9 ,  degrees 

o 

■ 

o 

0 . 2-. 

0.4 

0.6 

15 

1.41 

1.82 

2.64 

3.47 

22.5 

1.33 

1.66 

2.32 

2.98 

30 

1.25 

1.50 

2.00 

2.50 

45 

1 .12 

1.25 

1.50 

1.75 

60 

1.05 

1.10 

1.20 

1.30 

It  is  seen  that  the  amplification  of  wave  height  on  the  banks  is 
sensitive  to  both  bank  slope  and  center-channel  wave  amplitude. 
The  amplification  rises  sharply  with  increasing  amplitude. 

It  is  of  interest  that  Peregrine  (1968)  notes  the  triangular 
solution  in  a  particularly  elegant  form,  although  he  does  not  go 
on  to  provide  the  development  for  amplitude  variation  given  here. 
In  his  expression  for  ^  ,  he  takes  the  origin  at  the  channel 
vertex  for  which 


^(y.z)  -  (y2  +  z2)  /  4 

This  shows  clearly  that  the  isolines  of  ^  are  segments  of 
circles.  Since  in  a  triangle  the  banks  are  simply  radii  of  these 
circles,  the  solid  wall  boundary  condition  is  satisfied  by 
inspection.  Surprisingly,  all  cross-sectional  motion  must  be 
radial  from  the  channel  invert.  This  also  helps  to  make  clear 
why  all  triangular  sections  have  the  same  form  for  . 


No  other  analytical  solutions  have  been  found. 


r 


3.1.3  Numerical  Solutions 


We  have  implemented  a  numerical  model  based  on  the  foregoing 
theoretical  development  for  the  case  of  trapezoidal  channels. 
The  model  at  this  time  is  limited  to  determination  of  the  wave 
crest  profile  across  the  channel  to  investigate  bank  amplifica¬ 
tion  effects,  as  described  above  for  the  known  analytical 
solutions . 


It  is  first  necessary  to  solve  numerically  for  U*  throughout  the 
domain,  satisfying  the  Poisson  equation  and  boundary  conditions 
of  Equations  16.  Then  we  apply  Equation  29  to  determine  the 
surface  transverse  profile. 

We  have  adopted  a  simple  Simultaneous  Over-Relaxation  scheme, 
based  in  part  on  the  methods  described  in  Press,  et  al  (  1986). 
Since  we  are  restricting  attention  to  trapezoids  in  this  Phase  I 
numerical  effort,  we  need  model  only  half  the  channel,  from  the 
center  to  either  bank.  Then,  by  symmetry,  the  appropriate 
boundary  condition  on  y  ■  0  is  the.  same  ajs  for  a  solid  wall 
running  down  the  center  of  the  channel: 


at  y 


(30) 


Letting  j  and  k  denote  indices  in  the  y-  and  z-directions , 
respectively,  the  Gauss-Seidel  difference  expression  for  an 
interior  point  is  simply 

*  ,  n  n+1  n  n+1  o 

'J'j.k  -  I/*  *j  +  l,k  +  *j-l,k  +  'fj.k+l  +  Yj,k-1  “  Ax^  (31) 

in  which  we  assume  a  square  mesh.  Letting  the  over-relaxation 
weighting  factor  be  u,  t^e  updated  quantity  is 


i , *  / •  N  n 


(32) 


It  is  not  clear  what  the  optimum  choice  of  u>  might  be  for  our 
geometry.  But  in  view  of  the  speed  and  economy  of  desktop  PC's, 
expenditure  of  effort  to  determine  an  optimum  value  is  not 
worthwhile.  Instead,  we  have  simply  adopted  values  appropriate 
for  the  bounding  rectangular  domain,  assuming  these  to  be  a 
reasonable  guess  for  our  somewhat  truncated  region. 

Therefore,  we  take 

w  -  2  /  (1  +  /(I  -  P j) ]  (33) 

in  which  Pj  is  the  Jacobian  spectral  radius  given  by 

pj  m  1/2  [cos(it/J)  +  cos(tt/K)]  (34) 

with  J  and  K  being  the  (integer)  dimensions  of  the  gridr  Again, 
this  expression  assumes  a  square  mesh. 
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The  domain  is  sketched  in  Figure  10  a.  On  the  left  boundary,  we 
set  equal  to  the  old  values  one  column  to  the  right.  On  the 
bottom  boundary,  similarly,  we  set  equal  to  the  values  just 
above.  For  most  of  the  top  boundary,  we  set  ^  equal  to  the 
values  just  below  plus  a  constant  equal  to  the  space  step  times 
Co- 

On  the  right  (sloping)  boundary,  we  construct  a  perpendicular 
from  the  boundary  cutting  either  a  vertical  or  horizontal  side  of 
an  adjacent  cell  (depending  upon  bank  slope).  The  value  assigned 
to  the  boundary  point  is  the  value  interpolated  from  the  end¬ 
points  of  the  intercepted  side. 

Difficulty  comes  for  the  free-surface  condition  in  the  vicinity 
of  the  bank  for  gentle  bank  slopes  (small  9  ).  In  this  case, 
there  may  be  several  surface  points  for  which  there  does  not 
exist  a  valid  adjacent  point  on  the  second  row  (points  on  the 
second  row  lying  below  the  slope,  outside  the  domain)  on  which  to 
construct  the  surface  gradient.  This  case  is  sketched  in 
Figure  10b . 

After  some  experimentation,  we  have  adopted  the  following 
artifice  which  has  appeared  to  work  adequately.  From  the  surface 
point  of  interest,  point  A  for  example,  di-op  vertically  to  the 
bank  (point  B) .  Then  from  this  point,  construct  a  perpendicular 
to  the  bank  intersecting  the  free  surface  at  point  C.  By  the 
zero  gradient  condition  or  the  bank,  approximate  at  B  by  the 
interpolated  value  of  at  C  using  the  most  recent  values  on  the 
first  row>2  Then,  by  the  surface  gradient  condition,  add  the 
quantity  c0<5z  in  which  6z  is  the  distance  between  points  A  and  B. 

Finally,  at  the  conclusion  of  each  iteration,  the  value  of  ^(0,0) 
is  subtracted  from  every  grid  point,  thereby  fixing  the  arbitrary 
constant  in  ^  appropriately.  A  convergence  test  is  made  on  the 
trace  of  the  matrix,  with  termination  upon  achieving  a  relative 
change  of  less  than  10“^  (with  a  minimum  of  100  iterations, 
however ) . 

With  these  conditions,  we  have  successfully  simulated  a  range  of 
trapezoid  geometries  and  wave  heights.  The  source  code  for  the 
model  is  included  as  Appendix  A  of  this  report  (written  in  Pascal 
suitable  for  PC  applications).  Results  obtained  with  the  model 
will  be  shown  in  Section  4. 
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Figure  10a:  Numerical  Mesh 


C  A 


Figure  10b:  Sketch  Illustrating  Approximate  Surface 
Boundary  Condition  (see  text). 


3.2  Damping 

3.2.1  Frictional  Damping 


Frictional  damping  is  not  usually  much  of  a  consideration  for 
explosion  wave  propagation.  In  open  water,  radial  spreading 
causes  a  rapid  drop  in  height;  this  loss  is  compounded  by 
frequency  dispersion  if  the  depth  is  not  too  small.  Breaking  may 
also  be  encountered.  However,  in  the  case  of  the  surf-zone 
effect  (the  so-called  'Van  Dorn  effect'),  LeMehaute  (1-985)  has 
shown  that  friction  may  play  an  important  role  owing  to  the  long 
distances  of  propagation  over  a  gentle  slope  and  the  reduction  of 
radial  spreading  by  refraction.  Whether  shoaling  is  adequate  to 
maintain  a  destructive  breaker  or  whether  friction  may  be  just 
large  enough  to  keep  the  waves  sub-breaking,  could  be  a  critical 
question . 

For  waves  in  channels  it  is  possible  that  friction  may  be  of 
moderate  importance  since  radial  spreading  is  eliminated,  depths 
may  be  small,  and  breaking  is  restricted  to  limited  zones  at  the 
banks.  Of  course,  depending  upon  the  geometry,  other  decay 
mechanisms  such  as  channel  branching  -  and  irregularity  are  likely 
to  be  much  more  important  than  friction  in  causing  the  eventual 
demise  of  a  wave.  Overall,  friction  is  probably  a  secondary 
consideration,  although  still  of  interest. 

•  "  v 

The  long  waves  of  interest  here  may  be  approximately  cnoidal  in 
form,  or,  in  the  limit,  solitary.  This  indicates  that  the  rate 
of  loss  under  crests  will  be  much  greater  than  under  troughs 
owing  to  the  concentration  of  velocity  near  the  crests  and  to  the 
fact  that  the  dependence  on  velocity  is  quadratic. 

Consequently,  it  is  of  interest  to  assess  the  magnitude  of 
frictional  damping  from  the  standpoint  of  cnoidal  or  solitary 
waves  rather  than  sinusoidal  waves.  In  this  Phase  I  effort  we 
have  gone  part  of  the  way  toward  that  goal,  while  avoiding 
unnecessary  complexity,  by  starting  from  the  'hyperbolic  wave 
theory'  of  Iwagaki  (1967). 

Hyperbolic  waves  are  simply  a  practical  approximation  to  cnoidal 
waves  of  sufficient  length.  Cnoidal  wave  theory  is  expressed  in 
terms  of  the  elliptic  integrals  E(k)  and  K(k),  and  the  Jacobian 
elliptic  function  cn(0;k).  For  waves  with  sufficient  amplitude 
and  length,  the  parameter  k  is  very  nearly  unity  and  E(k)  is 
negligibly  greater  than  unity.  K(k),  however,  remains  finite, 
becoming  large  (the  solitary  wave  limit)  only  for  k  extremely 
near  unity;  departures  as  small  as  10“^0  are  not  insignificant. 

Iwagaki  (1967)  noted  that  a  useful  approximation  would  result  if 
we  simply  set  k  and  E(k)  to  unity,  and  determine  K(k)  as  a 
function  of  water  depth,  wave  height,  and  wave  period.  The 
Jacobian  elliptic  functions  reduce  to  their  hyperbolic-function 
limits  (hence,  the  name  'hyperbolic  waves'),  greatly  simplifying 
engineering  computations  at  little  sacrifice  of  accuracy.  In 
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particular,  the  function  cn(9;k)  becomes  sech(6)  and  sech(6)  is 
not  periodic.  This  poses  no  real  problem.  We  simply  imagine  a 
succession  of  waves  patched  together  where  the  troughs  should  be. 
At  these  patches,  the  wave  description  is  not  at  all  accurate 
(the  direction  of  flow  is  wrong,  for  example)  but  since  effects 
are  concentrated  at  the  crests  we  proceed  on  the  hope  that  the 
net  error  will  be  small.  In  this  Phase  I  effort  w§  further 
simplify  matters  (greatly)  by  truncating  all  expressions  for 
hyperbolic  wave  features  (particle  velocities,  energy  density, 
and  so  forth)  at  the  leading,  dominant  terms,  neglecting  in 
particular  powers  of  H/h. 

Given  these  assumptions,  we  proceed  as  follows.  Let  E  denote  the 
wave  energy  density  per  unit  area  averaged  over  the  length  of  the 
wave.  Assume  that  the  rate  of  energy  dissipation  per  unit  area 
is  given  by 

Et  -  E  -  kfPu3b  (35) 

Averaging  this  over  the  wavelength  we  write  the  energy  conserva¬ 
tion  equation  as  - 

_  • 

( EVg ) x  +  E  -  0  (36) 

in  which  Vg  is  the  group  velocity  and  ub  is^the  fluid  velocity  at 
the  bed;  kf  is  a  friction  coefficient  equal  to  f/8  in  which  f  is 
the  Darcy-Weisbach  friction  coefficient. 

From  hyperbolic  theory  we  use  the  approximation 

ub  -  /( gh )  H/h  sech2  ( 2Kx/L) 

in  which  K  is  approximated  by 

K  -  T/(3gH)/4h  (37) 

The  reader  will  note  that  this  is  essentially  the  solitary  wave 
approximation,  again  neglecting  H/h.  The  advantage  is  just  that 
we  retain  L  explicitly  for  purposes  of  defining  wavelength 
averages  (the  average  energy  density  over  a  solitary  wave  being 
zero ) . 

Substituting  into 

E  -  1/L  jL  E  dx 


we  find 

E  «  kfPH3  (g/h)3/tanhK  (8/3  +  4/3  sech2K  +  sech4K)/5K 
On  the  other  hand,  the  total  energy  density  is 

E  -  2/3  PgH2  (1-3  /(2K)]/K  r 
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from  which  we  can  evaluate  Ex  as 

Ex  -  PgHx/ A2  (A/H  -  1) 
in  which  A  is  the  constant  T/(3g)/4h. 

Substituting  these  expressions  for  Ex  and  Et  into  the  energy 
conservation  equation  gives  an  expression  for  Hx  which  can  be 
integrated  analytically  to  give  an  expression  relating  Hj,  H2. 
and  X  (the  subscripts  denoting  wave  height  at  the  beginning  and 
end  of  an  interval  X).  This  expression  is  a  cubic  in  /H2  which 
can  be  solved  in  the  form 

®2  +  P®  2  +  q  ■  0 

e2  -  /h2 

P  -  -  K2ei/(x6i  +  K20i  -  Hi) 

01  -  /H  x 

Hi  -  20h3/[kfT/(3g)] 

K2  -  15  h2  /2kf 

3  V>v 

q  -  Klei/(x6i  +  r20i  -  Ki) 

Graphs  showing  attenuation  computed  in  this  manner  are  shown  in 
Figures  11  -  13  for  initial  wave  amplitudes  (dimensionless)  of 
0.2,  0.4,  and  0.6,  respectively.  In  each  case,  H2  is  shown 

versus  reduced  propagation  distance  (X/L)  with  relative  wave 
length  as  a  parameter.  The  Darcy-Weisbach  coefficient  was  taken 
as  0.02  giving  a  value  for  kf  of  0.0025.  This  corresponds  to  a 
Manning  coefficient  of  about  0.028  in  a  water  depth  of  100  feet. 

Despite  the  appearance  of  these  curves,  it  will  be  seen  that  in 
practice  we  will  generally  encounter  dimensionless  propagation 
distances  on  the  order  of  10  to  20  so  that  typical  height 
reductions  will  be  on  the  order  of  only  10Z. 

3.2.2  Breaking  Dissipation 

We  have  found  no  information  directly  pertinent  to  defining  the 
breaking  zone  on  the  banks  and  the  associated  energy  losses. 
This  remains  an  important  item  for  future  research  including 
several  lines  of  approach.  For  example,  the  three-dimensional 
solitary  wave  theory  should  be  pursued  for  particle  velocities 
and  accelerations  in  order  to  define  a  breaking  criterion 
analogous  to  the  familiar  H/h  -  0.78.  This  condition  is  obtained 
from  the  McCowan  solitary  wave  theory  under  the  kinematic 
assumption  that  particle  velocity  at  the  crest  equals  wave 
celerity. 

Phenomenological  models  of  wave  breaking  should  also  be  inves- 
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Wave  Height  Reduction  (H2/H1) 


Figure  11:  Frictional  Height  Reduction  for  Initial 
Height  of  0.2. 
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Figure  12:  Frictional  Height  Reduction  for  Initial 
Height  of  0.4. 
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Figure  13:  Frictional  Height  Reduction  for  Initial 
Height  of  0.6. 


tigated.  The  non-saturated  breaker  theory  of  LeMehaute  (1962)  is 
an  example  of  such  a  theory  which  attempts  to  relate  breaking 
energy  dissipation  to  simpler,  underlying  concepts  such  as  rate 
of  energy  loss  in  an  analogous  hydraulic  jump. 

It  is  important,  however,  that  at  least  a  limited  experimental 
series  be  undertaken  to  give  empirical  guidance  in  this  area. 
Without  this,  any  theoretical  kinematic  or  dynamic  breaking 
criterion  must  be  considered  speculative. 

3.3  Channel  Branching  Effects 

The  final  topic  of  general  interest  is  the  height  loss  associated 
with  splitting  of  the  channel  into  branches.  This  can  be  an 
extremely  complex  problem  --  every  site  being  unique  --  if  the 
channel  junction  is  not  smaller  than  the  length  of  the  wave.  For 
short  waves,  the  methods  of  coastal  engineering  must  be  invoked, 
following  the  waves  by  ray-tracing  or  some  other  method  account¬ 
ing  for  refraction,  diffraction,  reflection,  and  so  forth.  No 
simple  general  rules  would  suffice. 

But  in  the  event  that  the  waves  are  reasonably  long  with  respect 
to  the  size  of  the  junction,  it  is  possible  to  make  a  very  great 
simplification  in  the  analysis.  In  this  approximation  --  the 
compact  junction  assumption  --  we  introduce  the  notions  of 
channel  admittance  and  impedance  (here  following  Lighthill, 
1978). 

If  the  junction  is  compact,  then  it  is  reasonable  to  assume  that 
each  branch  is  subject  to  much  the  same  surface  elevation,  or 
pressure  excess.  The  admittance  of  a  channel  is  defined  as  the 
ratio  of  volume  flow  to  pressure  excess  and  is  named  by  analogy 
with  the  electrical  admittance  function  which  is  the  ratio  of 
current  to  voltage  difference.  In  the  case  of  a  channel  of 
arbitrary  cross-section,  Lighthill  shows  that  the  admittance,  Y, 
is  given  by 


Y  -  A  /  (Pc) 

in  which  A  is  the  channel  cross-sectional  area,  P  is  density,  and 
c  is  wave  celerity.  The  impedance,  Z,  is  defined  as  the  recipro¬ 
cal  of  the  admittance. 

Now  consider  a  multiple  junction  --  assuming  continuity  of 
pressure  and  volume  flux  --  with  an  incident  wave,  f,  in  one 
branch.  Lighthill  shows  that  the  ratio  of  the  transmitted  wave 
height,  h  (which  is  the  same  in  all  secondary  branches),  to  the 
incident  height  is 


h/f  -  2Y !  /  ( Y i  +  EYi ) 


in  which  the  summations  are  taken  over  the  secondary  channels. 
Similarly,  the  reflected  wave  amplitude  ratio,  g/f,  is  given  by 


g/f  «  ( Y i  -  ZYi)  /  (YX  +  ZYi) 


In  these  expressions  we  assume  that  we  are  dealing  with  a  single 
wave  and  we  ignore  the  possibility  that  with  multiple  waves  the 
waves  reflected  from  one  branch  would  arrive  as  incidents  waves  at 
another  branch  so  that  superposition  would  have  to  be  considered. 

Lighthill  treats  the  case  of  periodic  waves  through  multiple 
branches  by  a  straight-forward  extension  of  the  above  ideas. 

In  our  case,  however,  we  will  be  content  to  consider  each  wave  as 
a  separate  entity. 

Despite  the  approximate  nature  of  this  treatment,  it  is  felt 
adequate  for  our  purposes.  Its  primary  limitation  is  the 
assumption  that  junctions  are  compact.  But  even  in  the  case  that 
this  is  not  strictly  true,  the  above  expressions  still  give  a 
useful  first  order  estimate  which  can  only  be  improved  with 
j  considerable  effort. 
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4.  PREDICTIONS  AND  COMPARISONS  WITH  OBSERVATIONS 
4.1.  Parametric  Presentation  of  Numerical  Results 

We  have  exercised  the  numerical  model  described  in  the  previous 
chapter  for  a  range  of  trapezoidal  channel  sizes  and  shapes,  and 
for  a  range  of  center  wave  amplitudes.  s 

The  basic  results  are  shown  in  the  following  five  four-part 
figures,  Figures  14-18,  showing  results  for  bank  slopes  of  15, 
22.5,  30,  45,  and  60  degrees,  respectively.  In  each  figure  we 
show  four  curves,  corresponding  to  relative  wave  amplitudes  at 
the  centers  of  the  channels  of  0.1,  0.2,  0.4,  and  0.6.  The  four 
parts  of  each  figure,  designated  a-d,  present  data  for  bottom 
half-widths  of  1.0,  1.5,  2.0,  and  5.0,  respectively. 

These  figures  clearly  show  the  expected  trends.  Amplification 
increases  with  increasing  channel  width,  with  increasing  central 
wave  height,  and  with  decreasing  bank  slope.  Amplifications  of  2 
to  3  appear  typical. 

In  order  to  demonstrate  the  dependence  on  bottom  width,  holding 
other  factors  constant,  we  present  Figure  19.  This  figure  shows 
the  amplification  at  the  shoreline  only,  as  a  function  of  bottom 
width,  holding  bank  slope  constant  at  15  degrees.  The  result  is 
clearly  that  within  the  range  of  widths  cbhsidered  here,  ampli¬ 
fication  is  not  sensitive  to  width. 

This  finding  is  very  comforting  in  view  of  the  limitation  of  the 
theory  to  channels  which  are  not  too  wide.  The  indication  is 
that  for  widths  as  great  as  five,  or  so,  spurious  amplification 
associated  with  excessive  width  is  not  a  problem. 

In  Figure  20,  we  fix  the  bottom  width  at  1.5  and  show  amplifica¬ 
tion  at  the  banks  versus  bank  slope.  Here  we  see  a  strong 
dependence  on  slope,  as  expected.  We  also  see  a  deficiency  of 
the  free-surface  boundary  condition  treatment  in  the  case  of  a 
very  gentle  slope.  As  already  discussed,  in  our  simple  square 
mesh  scheme,  several  surface  points  are  estimated  using  an 
approximate  technique  when  the  slope  is  15  degrees.  While  that 
technique  gives  a  reasonable  estimate  and  does  not  destroy 
convergence,  it  does  appear  to  underestimate  the  bank 
amplification  by  about  10%  in  this  case. 

In  order  to  see  this,  we  have  plotted  the  values  corresponding  to 
the  analytical  solution  for  a  15  degree  triangle;  these  points 
give  the  dashed  extensions  of  the  curves  to  15  degree  slope. 
Noting  from  Figure  19  that  bottom  width  is  not  a  sensitive 
parameter,  it  is  reasonable  think  that  this  analytical  solution 
will  be  close  to  the  desired  trapezoidal  value  for  this  width. 

Indeed,  Figure  20  shows  that  the  points  corresponding  to  the 
triangle  solution  do  appear  to  fall  at  about  the  expectea  values. 
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Figure  14b:  Amplitude  Variation  for  Indicated  Condition. 
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Figure  14c:  Amplitude  Variation  for  Indicated  Condition. 


r 


49 


I  Hm  am  • 


Figure  14d:  Amplitude  Variation  for  Indicated  Condition. 
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Figure  15b:  Amplitude  Variation  for  Indicated  Condition 


Figure  15c:  Amplitude  Variation  for  Indicated  Condition. 
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Figure  15d:  Amplitude  Variation  for  Indicated  Condition. 
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Figure  16a:  Amplitude  Variation  for  Indicated  Condition. 


Figure  16b:  Amplitude  Variation  for  Indicated  Condition. 
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Figure  16c:  Amplitude  Variation  for  Indicated  Condition. 
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Figure  16d:  Amplitude  Variation  for  Indicated  Condition. 
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Figure  17a:  Amplitude  Variation  for  Indicated  Condition. 
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Figure  17d:  Amplitude  Variation  for  Indicated  Condition. 
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Figure  18a: 


Amplitude  Variation  for  Indicated  Condition 
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Figure  18b:  Amplitude  Variation  for  Indicated  Condition. 
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Figure  18d:  Amplitude  Variation  for  Indicated  Condition. 
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Figure  19:  Variation  of  Amplitude  at  Bank  with  Bottom  Width 


Amplification 


Figure  20:  Variation  of  Amplitude  at  Bank  with  Bank  Slope. 
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4.2.  Comparisons  with  Observations 
4.2.1.  Data  of  Peregrine 


We  have  already  discussed  Peregrine's  (1969)  data  in  connection 
with  Figure  8.  Suffice  it  to  say  here  that  we  have  carefully 
scaled  the  surface  profile?  of  his  two  figures  to  determine  his 
observed  bank  amplif  icati-ons ,  and  have  then  compared  the  results 
of  our  numerical  model  with  those  values.  In  every  case,  our 
model  agrees  as  well  with  the  data  or  better  than  does  Pere¬ 
grine's  own  computed  points. 


4.2.2.  Lituya  Bay  Landslide  Wave 

Comparisons  here  are  really  speculative  but  nevertheless  of 
interest.  While  we  have  some  idea  of  the  peak  water  levels  from 
Figure  1,  we  do  not  have  observations  of  the  waves  at  the  center 
of  the  Bay.  However,  we  can  borrow,  the  numerical  estimates  of 
Chiang,  et  al  (1981)  who  simulated  this  case. 

The  basic  results  of  Chiang,  et  al,  are  summarized  in  the 
following  table  where  the  nine  station  locations  are  as  indicated 


in  Figure  la. 

Station  # 

Model  Peak  HeiRht 

Nearby  Observed 

1 

200  ft 

>  1000 

2 

120 

350 

3 

80 

130 

4 

60 

100 

5 

30 

100 

6 

50 

90 

7 

20 

— 

8 

40 

85 

9 

25 

30 

Clearly,  at  Stations  1  and  4,  the  comparisons  are  not  appropriate 
since  what  is  actually  occurring  is  a  normal  sort  of  runup  in 
those  areas.  Station  7  is  in  the  shadow  of  the  islan'd  and  is 
difficult  to  interpret.  But  for  all  other  points  where  fwe  might 
expect  some  sort  of  bank  amplification  we  do,  in  fact,  find  one. 
If  the  model  values  are  assumed  to  be  reasonable  estimates,  then 
the  apparent  amplification  is  about  2  to  3  typically. 
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It  does  appear  reasonable  to  think  the  computed  values  are  valid. 
To  be  greatly  off,  it  would  be  necessary  to  have  mis-es  t  i  ma  t  e  d 
the  slide  mass  by  a  considerable  amount.  Since  the  slide  energy 
goes  like  the  mass  and  since  wave  height  goes  like  tt*e  square 
root  of  energy,  a  100%  error  in  slide  mass  estimation  would 
correspond  to  about  a  40%  error  in  computed  heights.  It  is 
doubtful  that  the  indicated  values  are  that  much  in  error. 

Although  this  example  is  very  indirect  and  somewhat  speculative, 
it  is  nevertheless  welcome  in  view  of  the  general  lack  of  data  of 
any  other  sort.  It  does  appear  to  be  consistent  with  the  general 
idea  that  amplifications  on  the  order  of  a  factor  of  two  are 
typical. 


5.  EXAMPLES 


In  this  section  we  give  one  detailed  example  of  the  foregoing 
methods  in  estimating  explosion  wave  vulnerability,  and  one  more 
speculative  example. 
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Assume  the  incident  channel  has  a  depth  of  about  500  fe;et  and  a 
cross-sectional  area  of  about  3X10®  f  1 2 .  Let  the  left  branch 
have  a  depth  of  about  500  feet  also,  with  an  area  of  about  7X10'; 
let  the  right  branch  have  a  depth  of  about  300  feet  and  an  area 
of  about  2X106. 


A 


The  second  example  is  a  rough  illustration  only.  Figure  23  is 
reproduced  from  Jane’s  Defence  Weekly,  Sept  6,  1986,  and  shows  a 
satellite  photo  of  a  supposed  Soviet  submarine  base.  According 
to  the  accompanying  article,  the  facility  is  for  'Typhoon'  class 
submarines.  The  article  further  suggests  that  the  quays  visible 
in  the  photo  (within  the  box)  are  on  the  order  of  200  meters  in 
length . 

Taking  this  as  a  rough  guide,  we  estimate  the  total  width  of  the 
channelized  geometry  to  be  about  2400  feet.  Assuming  a  slope  of 
about  15  degrees  and  a  depth  of,  say,  250  feet,  we  arrive  at  an 
estimate  of  bottom  half-width  (dimensionless)  of  about  1.5. 


t 

c 
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Figure  23:  Aerial  Photo  of  a  Soviet  Submarine  Base  in  a 
Channelized  Geometry  (Kola  Peninsula;  Photo 
Reproduced  from  Jane's  Defence  Weekly, 

Sept.  6,  1986). 
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6.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  Phase  I  effort  has  summarized  the  major  known  features  of 
large  amplitude  cnoi da  1  or  solitary  waves  in  channels,  and  has 
provided  simple  formulas  to  estimate  wave  generation  by  both 
explosions  and  landslides.  The  salient  points  are: 

In  a  uniform  channel,  waves  may  propagate  long  distances 
with  little  loss  owing  to  the  smallness  of  the  friction 
effect  and  the  confinement  afforded  by  a  channel. 

Opposing  influences  include  wave  diminution  by  channel 
branching,  channel  irregularity,  and  breaking  in  the  shallow 
wedge-like  segment  of  the  crest  propagating  along  the  dry 
bank . 

In  the  case  of  sufficiently  long  waves,  channel  branching 
effects  are  reasonably  easily  estimated,  and  channel 
irregularities  become  secondary. 

Breaking  losses  are  still  unknown.  In  wide  channels  with 
moderate  slopes  the  losses  will  be  small  since  the  breaking 
zone  is  expected  to  be  very  limited.  In  any  case,  the 
breaking  zone  coincides  with  the  target  zone  so  that  this 
form  of  dissipation  is  not,  in  a  sense,  'wasted.' 

There  may  be  significant  enhancement  of  wave  height  on  the 
banks,  a  sort  of  'lateral  shoaling'  in  the  shallower  regions  ; 

on  the  slope  (although  the  effect  extends  even  beyond  the  r 

bank  into  the  flat-bottom  region  of  a  trapezoidal  channel).  | 

Figure  24  shows  the  approximate  upper  limit  of  the  amplifi-  = 

cation  factor  as  a  function  of  bank  slope  and  center  channel  c 

amplitude.  The  more  extreme  portions  of  this  figure  would 
be  greatly  modified  by  breaking  effects. 

Additional  effects,  not  specifically  discussed  earlier,  may  also 
play  a  significant  role  in  the  wave  description.  For  example,  a 
channel  with  gradually  changing  depth  or  width  may  amplify  (or 
attenuate)  the  wave  owing  to  familiar  mechanisms  of  shoaling  and 
compression  of  wave  orthogonals  (coincident  with  the  banks).  In 
particular,  there  may  exist  important  sites  for  which  wave 
generation  could  occur  in  relatively  deeper  water  (say,  offshore 
or  in  an  embayment)  and  then  propagate  into  the  channel 
(upstream).  In  this  way,  the  center  channel  wave  amplitude  might 
be  much  greater  than  could  be  achieved  by  a  detonation  within  the 
channel  itself.  This  would  have  the  double  effect  of  increasing 
the  bank  amplification  factor  as  well  (moving  up  to  a  higher 
curve  in  Figure  24).  . 

This  Phase  I  effort  was  primarily  intended  to  be  ar  problem 
review  and  discussion.  In  Phase  II  we  recommend  that  a  research 


78 


alt  U.S.  Uove 


1 


2 


10 


11 


12 


3  4  5  6  7  8  9 

Reciprocal  of  bank  slope  (cote) 


Figure  24:  Approximate  limit  of  bank  amplification  as  a  function  of  bank  slope. 

These  curves  represent  the  triangular  cross-section  analytical 
solution,  neglecting  the  additional  enhancement  obtained  in 
wider  channels.  r 
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effort  be  undertaken  to  extend  our  understanding  of  these 
questions.  In  particular,  we  suggest  the  following. 

The  three-dimensional  wave  theory  should  be  further  pursued 
to  eliminate  the  limitations  discussed  earlier.  That  is,  it 
should  be  extended  to  handle  both  wide  channels  (single 
banks  in  the  limit)  and  gentle  slopes. 

The  numerical  implementation  given  here  should  be'  extended 
to  handle  irregular  channel  geometry  and  to  improve  the 
treatment  of  the  boundary  conditions. 

The  numerical  model  should  be  exercised  for  a  variety  of 
channels  of  interest. 

The  three-dimensional  theory  should  be  pursued  to  derive 
expressions  for  fluid  particle  velocities  and  accelerations. 

Using  these  results,  breaking  criteria  should  be  sought. 
For  example,  particle  velocities  and  accelerations  along  the 
portion  of  the  crest  lying  on  the  dry  bank  should  reveal  the 
occurrence  of  a  critical  condition  indicating  incipient 
disintegration  of  the  wave  (either  kinematic  or  dynamic). 

A  supplementary,  small  experiment  s'!  program  should  be 
performed  to  guide  and  verify  analytical  work,  especially  in 
the  matter  of  breaking  criteria. 

Based  on  all  the  above,  improved  methods  of  estimating  bank 
heights  should  be  developed. 

Similarly,  expressions  should  be  derived  relating  gross  wave 
and  channel  characteristics  to  over-bank  dynamic  effects 
(wave  forces  on  bank  structures). 

The  implications  of  these  results  should  be  assessed  for 
wave  behavior  in  non-channelized  geometries.  For  example,  a 
so  1  i  t a r y -1 i ke  waveform  in  a  complex  three-dimensional 
geometry  may  at  some  points  interact  with  slopes  in  a  manner 
giving  rise  to  amplification.  This  might  be  ultimately 
related  to  the  problem  of  behavior  in  a  Mach-stem  zone. 

The  developed  methods  should  be  implemented  to  assess  the 
vulnerability  of  selected  sites  (worldwide)  to  both  ex¬ 
plosion  and  landslide  effects  using  the  most  appropriate 
generation  descriptions  available. 


r 


80 


/ 


J 


I 


7.  REFERENCES 


Benet,  F.  and  J.A.  Cunge,  1971,  "Analysis  :>  f  Experiments  on 
Secondary  Undulations  Caused  by  Surge  Waves  in  Trapezoidal 
Channels,"  Journal  of  Hydraulic  Research,  v9. 

Benjamin,  T.B.  and  M.J.  Lighthill,  1954,  "On  Cnoidal  Waves  and 
Bores,"  Proceedings  of  the  Royal  Society,  A224, 

Chiang,  W.L.,  Divoky,  D.,  Parnicky,  P.,  and  W.  Wier  (1981), 
"Numerical  Model  of  Landslide-generated  Waves,"  Bureau  of 
Reclamation  0-07-8 1 -SO  1 34 . 

Divoky.  D.  and  John  Lane,  1973,  "Edge  Bores,"  Office  of  Naval 
Researc.  '00014-72-C-0366 ,  Tetra  Tech  Rept  TT-P-282-73- 1 . 

Fenton,  J.D.,  1973,  "Cnoidal  Waves  and  Bores  in  Uniform  Channels 

of  Arbitrary  Cross-Section,"  Journal  of  Fluid  Mechanics,  v58, 
part  3,  pp  417-434. 

Iida,  K.,  D.  Cox,  and  G.  Pararas-Carayannfs,  1967,  "Preliminary 
Catalog  of  Tsunamis  Occurring  in  the  Pacific  Ocean,"  Hawaii 
Institute  of  Geophysics,  HIG-67-10. 

Iwagaki,  Y.,  1968,  "Hyperbolic  Waves  end  their  Shoaling," 

Proceedings,  11th  Conference  on  Coastal  Engineering. 

Korteweg,  D.J.  and  G.  de  Vries,  1895,  "On  the  Change  of  Form  of 
Long  Waves  Advancing  in  a  Rectangular  Canal  and  on  a  New  Type  of 
Long  Stationary  Wave,"  Philosophical  Magazine,  v39,  n5. 

Kranzer,  H.  and  J.  Keller,  1959,  "Water  Waves  Produced  by 
Explosions,"  Journal  of  Applied  Physics,  v30,  n3. 

LeMehaute,  B.,  1962,  "On  the  Nonsaturated  Breaker  Theory  and  the 

Wave  Runup,"  Proceedings,  8th  Conference  on  Coastal  Engineering. 

LeMehaute,  B.,  Divoky,  D.,  and  A.  Lin  1968,  "Internal  Charac¬ 
teristics  of  Explosion-generated  Waves  on  the  Continental  Shelf," 
Defense  Atomic  Support  Agency,  DASA  01 -67-C-0099 . 

LeMehaute,  B.,  1971,  "Theory  of  Explosion-generated  Water  Waves," 
in  Advances  in  Hydroscience,  v7.  Academic  Press,  N.Y. 

LeMehaute,  B.,  1985,  "Explosion  Generated  Waves  and  Submarine 

Response,  Volume  1,"  DNA,  Contract  DNA-001-84-0053 ,  February. 

Lighthill,  J.,  1978,  Waves  in  Fluids,  Cambridge  University  Press, 
Cambridge . 

Miller,  D.J.,  1960,  "Giant  Waves  in  Lituya  Bay,  Alaska,"  USGS 

Professional  Paper  354-C. 


81 


Noda,  E.K.  ,  1969,  "Theory  of  Water  Waves  Generated  by  a  Time- 
Dependent  Boundary  Displacement,"  Technical  Rept.  HEL  16-5, 
University  of  Calif.,  Berkeley. 

Peregrine,  D.H.,  1968,  "Long  Waves  in  a  Uniform  Channel  of 
Arbitrary  Cross-Section,"  Journal  of  Fluid  Mechanics,  v32,  part 
2,  pp  353-365. 

Peregrine,  D.H.,  1969,  "Solitary  Waves  in  Trapezoidal  Channels," 
Journal  of  Fluid  Mechanics,  v35,  part  1,  pp  1-6. 

Peters,  A.S.,  1966,  "Rotational  and  Irrotational  Solitary  Waves 
in  a  Channel  with  Arbitrary  Cross-Section,"  Communications  on 
Pure  and  Applied  Mathematics,  vl9,  pp  445-471. 

Press,  W.H.,  B.P.  Flannery,  S.A.  Teukolsky,  and  W.T.  Vetterling, 
1986,  Numerical  Recipes,  The  Art  of  Scientific  Computing, 
Cambridge  University  Press,  Cambridge. 

Raney,  D.,  and  H.L.  Butler,  1975,  "A  Numerical  Model  for  Predict¬ 
ing  the  Effects  of  Land s 1 ide -Ge ne r a t ed  Water  Waves,"  US  Army 
Corps  of  Engineers,  Waterways  Experiment  Station,  Report  H-75-1. 

Slingerland,  R.,  and  B.  Voight,  1982,  ftTvaluating  Hazard  of 
Landslide-Induced  Water  Waves,"  ASCE  Journal  WW4 ,  November. 

Strange,  J.N.,  1955,  data  contained  in  "Effects  of  Explosions  in 
Shallow  Water,"  US  Army  Corps  of  Engineers  Waterways  Experiment 
Station,  Technical  Memorandum  2-406. 

Ursell,  F.  ,  1953,  "The  Long-Wave  Paradox  in  the  Theory  of  Gravit  y 
Waves,"  Proceedings  of  the  Cambridge  Philosophical  Society,  v49. 

Wang,  S.,  R.  Wade,  and  W.  Wier,  1977,  "Vulnerability  of  Surfact 
Effect  Vehicles  to  EGWW,"  prepared  for  Office  of  Naval  Research, 
Tetra  Tech  Rept.  TC-645. 

Wiegel,  R.  1964,  Oceanographical  Engineering,  Prentice-Hall,  N.J. 


82 


J 


APPENDIX  A 


C 


this  copy  meoe  at  U.S.  Govei 


Listing  of  Channel  Wave  Program  CHANWAVE . PAS ,  page  1 

<*V->  {$R+ }  <$U->  {$C+>  {*K-> 

program  bankWaves; 

Procedure  Pause; 

YAT 

junk  :  char; 

Begin 

writer  ..  Press  Any  Key  To  Continue  ..  ’ : 45 ) ; 

read(kbd, junk) ; 

writeln; 

find; 


{=======================  TIMER  MODULE  ======================= 

const 

nTimers  =  5; 

YAX 

beginTime:  arrayTl.  .nTiaere-)  of  real; 
procedure  initialiseTimers; 

YAX 

i:  byte; 
const 

infinity  =  1.0ES0; 
begin 

far  i  ::  1  tfl  nTimera  dp  beginTime[i]  :=  infinity 
end:  v 


function  currentTime:  real; 
type 

regPack  =  rflfiprd 

ax,  bx,  cx,  dx,  bp,  ai,  di,  ds,  es,  flags:  integer; 

find; 

YAX 

regs:  regPack; 
begin 

with  regs  do 
begin 

ax  :=  S2C00; 
msdos(regs) ; 

currentTime  :=  lo(dx)/100  ♦  hi(dx)  ♦  60  *  (lo(cx)  +  60  *  hi(cx)) 
end 

find; 


PXftCfld-UXfi  beginTimer (n:  byte); 

bag  in 

if  n  in  [1.. nTimers]  than  beginTime[n]  :s  currentTime 

find; 


fnnpiinn  elapsedTime(n:byte) :  real; 
const 

infinity  =  1.0E30; 
begin 

if  n  in  [1.. nTimers]  then  elapsedTime  :=  currentTime  -  beginTime[n) 
else  elapsedTime  :=  infinity 

end ; 


{======================  end  0F  TIMER  MODULE  =====================) 


Listing  of  Channel  Wave  Program  CHANWAVE . PAS ,  page  2 


type 

anyStr  =  stringr2551 ; 
str2  r  string! 21 ; 
str3  =  aitingC3]; 

const 

omegaFactor  =  1.0; 
minlterations  =  100; 
maxlterations  =  2000; 
maxRows  =  21; 
maxCols  s  151; 
epsilon  =  1.0E-7; 
infinity  =  1.0E30; 
tolerance  =  1.0E-6; 

var 

response:  char; 

index,  row,  col,  nRows,  nCols,  iterCount, 

showSteps,  zSteps ,  baseCol,  midCol,  coll4,  col34,  triCol:  integer; 

midMax,  lastCol:  g£rgg[l • .maxRowsJ  sH  integer; 

depth,  width,  alpha,  oldSum,  newSum,  topWidth,  heightO, 

delY,  delZ,  beta,  betaSqr,  betaSqlnv,  spectralRadius ,  omega, 

aZero,  bZero,  delySqr,  delSurf,  factorl,  factor2, 

tanAlpha,  cot Alpha,  aurfGrad,  bZeroPrime,v uNew,  theta,  omegaOpt, 

coeffl,  coeff 2,  duration,  oldCorner,  newCorner,  convergeTest :  real; 

u:  nrcazt  1  •  -  maxRows ,  l..maxCols]  si  real; 

height:  .  .maxCols]  at  real; 

OK,  ahowResulta,  converged,  save:  boolean; 

onBoundary:  g££gy[l . .maxRows]  q1  boolean; 

junction  tan (theta:  real):  real; 
begin 

tan  :s  s in (theta) /cos (theta) 

find; 


Er0CfidUC.fi  bell(freq,  time:  integer); 
begin 

sound (freq) ; 
delay (time) ; 
noSound ; 

find; 


ECggfldUEfl  wait (message:  anyStr); 
YAC 

ch:  char; 
begin 

write(message) ; 
repeat  until  keypressed; 
read( kbd , ch) 
end: 


procedure  get Input; 


Listing  of  Channel  Wave  Program  CHANWAVE . PAS ,  page  3 

clrscr ; 

write (’Enter  depth  of  the  channel:  *:55); 
readln( depth) ; 

write (’Enter  number  of  steps  in  depth  direction:  ’:55); 
read In ( zSteps ) ; 

write(’Enter  bottom  width  of  channel:  ’:5S); 
readln(width) ; 

write( ’Enter  bank  inclination  (degrees  from  horizontal):  ’:55); 
readln( theta) ; 
alpha  : =  90  -  theta; 

If  (theta  <=  0)  ax.  (theta  >90)  than 
begin 

bell (500,  500); 

writeln( ’Bank  slope  out  of  legal  range. . .exiting.  .  .  ’ ) ; 
wait (’'... press  any  key ...’); 
halt 
end; 

write(’Enter  wave  amplitude  at  left  boundary:  ’:55); 
readln(heightO) ; 
showResulta  :  =  false; 

write ('... Want  to  see  intermediate  results...  Y/M  ?:  * : 55 ) ; 

repeat 

read(kbd,  response); 
response  : =  upcase( response) 
until  response  in  [ ’ Y ’ ,  ' N ’ ] ; 
writeln( response ) ; 
if  response  =  ’Y’  then 

begin  'v 

showResulta  : =  true; 

write( ’.. .After  every  how  many  iterations  ?:  ' : 55 ) ; 
read In ( showSteps ) ; 
end; 

write (’Save  results  to  a  text  file...  Y/N  ?:  *:55); 

roreat 

read ( kbd ,  response ) ; 
response  : =  upcase ( response ) 
until  response  in  C’Y’,  ’  N’]; 
write In (response) ; 
save  : =  response  =  ’Y’ ; 

and; 

procedure  initialise; 
var 

zRow,  widthRow,  zPtBelow,  yPtBelow,  zPtRight,  yPtRight:  real; 
begin 

converged  : =  false; 
iterCount  :=  0; 
oldSum  :s  0; 
newSum  : =  0 ; 
oldCorner  :=  0; 
newComer  :=  0; 
convergeTest  : =  1.0; 
delZ  :=  depth  /  zSteps; 
delY  :s  delZ; 
beta  :=  delY/delZ; 
delySqr  :=  delY  *  delY; 
nRows  : =  zSteps  ♦  1 ; 
alpha  :  =  alpha  *  pi  /  180; 
tanAlpha  :=  tan( alpha); 


t  A  A  on 


Listing  of  Channel  Wave  Program  CHANWAVE . PAS ,  page  4 

if  theta  =  90  than  cotAlpha  :  =  infinity  else  cotAlpha  :=  1  /  tanAlpha; 

factorl  :=  cotAlpha  /  beta; 

factor2  :=  beta  *  tanAlpha; 

topWidth  :=  width  +  depth  *  tanAlpha; 

nCols  :=  trunc( topWidth  /  delY  +  epsilon)  +  1; 

for  col  : =  1  to  nCols  dfi  height[col]  :  =  0; 

for  row  : =  1  to  nRows  da  for  col  : =  1  to  nCols  dfi  u[row,col]  :  =  0; 
for  row  : =  1  to  nRows  do 

zRow  :=  delZ  *  (nRows  -  row); 
widthRow  :=  width  +  zRow  *  tanAlpha; 

lastColfrow]  :=  trunc( widthRow  /  delY  +  epsilon)  +  1; 

onBoundary[row]  :=  ((widthRow  -  ( lastCol [row]  -  1)  *  delY)  <  epsilon); 

and: 

for  row  : =  1  to  nRows  do 

for  col  : -  1  to  lastCol[row]  do 

u[row,col]  :=  tanAlpha  *  sqrt( topWidth)  * 

(1  -  sqrt(l  -  (col  -  1)  *  delY  /  topWidth))  -  0.4  *  (row  -  1)  *  delZ; 
baseCol  :=  trunc( width  /  delY  +  epsilon); 
midCol  :s  baseCol  div  2; 
col 14  : =  baseCol  div  4; 
col34  : =  3  *  coll4; 
triCol  :  =  (baseCol  +  nCols)  div  2', 
for  row  : =  1  to  ( nRows  -  1 )  do 
begin 

col  :=  baseCol; 
repeat 

yPtBelow  :s  delY  *  (col  -  1); 
zPtBelow  :=  delZ  *  (nRows  -  row  -  1); 

OK  :=  ((width  +  zPtBelow  *  tanAlpha  -  yPtBelow)  >=  -epsilon); 
zPtRight  :s  delZ  *  (nRows  -  row); 
yPtRight  :s  delY  *  col; 

if  OK  then  OK  :=  ((width  +  zPtRight  *  tanAlpha  -  yPtRight)  >s  -epsilon); 
col  :=  col  *  1 
until  not  OK; 
midMaxCrow]  : =  col  -  2 

and: 

midMax [nRows]  :s  baseCol; 
writeln; 

writeln( nRows,  '  rows,  nCols,  ’  columns.’); 
writeln; 

for  row  ::  1  to  nRows  do 

writeln(row:5,  midMax [row] : 10,  lastCol[row] : 10,  onBoundary[row] : 10) ; 
writeln; 

wait ( ’ . . . press  any  key ...’); 
writeln;  writeln; 
betaSqr  : =  beta  *  beta ; 
betaSqlnv  : =  1  /  betaSqr ; 
aZero  : =  width  ♦  0.5  *  tanAlpha; 
bZero  :=  topWidth; 
surfQrad  :=  aZero  /  bZero; 
delSurf  :=  surfGrad  *  delZ; 
bZeroPrime  :=  tanAlpha; 
spectralRadius  : = 

(cos(pi/(nRows-l ) )  +  betaSqlnv  *  cos ( pi/( nCols- 1 )) )  /  (1  ♦  betaSqlnv); 
omegaOpt  :=  omegaFactor  *2/(1  ♦  sqrt(l  -  sqr( spectralRadius ))) ; 
omega  ■=  omegaOpt; 

writeln( ’ Columns  printed:  ',  ’l’:5,  coll4:5,  midCol:5, 
col34:5,  baseCol: 5,  triCol: 5,  ’  plus  points  on  slope.’); 
writeln; 
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coeffl  :=  omega  /  (2  *  (1  ♦  betaSqr)); 
coeff2  :=  1  -  omega; 
and: 


procedure  topRow; 

procedure  upperLef tCorner ; 
begin 

ufl.l]  :=  u[l , 2] 

gj od; 

procedure  midTop; 
begin 

for  col  :=  2  midMax[l]  dp  u[l,col]  :=  u[2,col]  +  delSurf 
and; 

procedure  upperRight Wedge; 

Yflr 

dalHaight,  factor:  real; 

function  useTop(col:  integer):  boolean; 
begin 

delHeight  :=  1  -  delY  *  (col  -  baaeCol  -  1)  *  cotAlpha; 
factor  : =  delHeight  *  cotAlpha  /  delY; 
uaeTop  :s  factor  <=  1; 

grid; 

begin 

£stX  col  :=  (midMax[l]  ♦  1 )  to  laatCol[l] 
if  useTop(col)  then 

u[l,colj  :=  u(l,col]  ♦  factor  *  (u[l, col-1]  -  u[l,col])  + 
surf Grad  *  delHeight 

elae 

ud.col]  :*  u[l, col-1]  ♦  (delHeight/delZ  -  factor2)  « 

(u[2, col-1]  -  u[l, col-1])  ♦  aurfOrad  *  delHeight; 
H  delHeight  <  epailon  then 

u[l,col]  :r  3  *  (u[l,col-l]  -  u[l, col-2])  ♦  u[l, col-3] 

end: 

begin 

upperLef tCorner ; 
midTop ; 

upperRightWedge 

gad; 


procedure  mldRowt row:  integer); 
procedure  leftPoint; 
begin 

u[row,l]  :s  u[row,2] 

gad; 

procedure  midPointa; 
begin 

for  col  : =  2  to  midMax[row]  do 
begin 

uHew  :=  u[roe,col+l]  ♦  u[row, col-1] 

♦  betaSqr  *  (u[row*l ,col]  ♦  u[row-l ,col] )  -  delySqr; 
u[ row, col]  : =  coeffl  *  uNew  *  coeff2  *  u[row,col]; 
and 

gad; 

procedure  rightPoints; 
haAlD 

lax  col  :=  (midMax[row]  ♦  1 )  la  lastCol[row]  ifi 
u[row,col]  :=  u[row, col-1]  ♦ 

(1  /  (1  ♦  cotAlpha  /  beta))  *  (u[row-l ,col]  -  u[row, col-1] ) 

end: 

begin 
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leftPoint ; 
midPoints; 
rightPoints 
And; 


procedure  bottomRow; 

procedure  lowerLef tCorner ; 
begin 

u[nRows,l]  :=  u[nRows,2] 
and: 

procedure  midBottom; 
begin 

for  col  : =  2  to  mi dMax[ nRows]  dfl  u[nRows,col]  :s  u [ nRows -1 ,col] 

find; 

procedure  lowerRightCorner; 

XAT 

uAs If Bottom,  uAsIfSide:  real; 
begin 

uAslf Bottom  :=  u[nRows-l , lastCol [nRows] ] ; 
if  theta  <  45  then 

uAsIfSide  :s  u[nRows-l,lastCol [nRows]]  + 
factorl  *  (u[nRows-l, lastCol [nRows] -1]  -  u[nRows-l , lastCol [nRows] ] ) 
else 

uAsIfSide  :s  u[ nRows, lastCol [nRows] -1]  + 
factor2  *  (u[nRows-l , lastCol [nRows ]-l]  -  u[ nRows, last Col [nRows] ]) ; 
u[ nRows, lastCol [nRows]]  :s  0.5  *  (uAs If Bottom  ♦  uAsIfSide) 
end; 

begin  'v 

lowerLef tCorner ; 
midBottom; 
lowerRightCorner 
find; 


procedure  normalise; 

XOX 

uZero:  real; 
begin 

uZero  :=  u[l,l]; 

oldSum  :s  newSum; 

newSum  : =  0 ; 

oldComer  :s  newCorner; 

iflX  row  : =  1  to  nRows  do 

for  col  : =  1  to  lastCol [row]  do 

begin 

u[row,col]  :=  u[row,col]  -  uZero; 
newSum  : =  newSum  ♦  abs(u[row,col] ) 

find; 

newCorner  : =  u( nRows, lastCol [nRows]] 
end; 


function  psiA:  real; 

XfiX 

sum:  real; 
row,  col:  integer; 
begin 

sum  : =  0 ; 

for  row  : =  1  to  ( nRows  -  1 )  do 
IflX  col  : =  1  to  (nCols  -  1 )  do 

sum  :=  sum  +  u[row, col]*u[row, col  +  1 ] +u[row+l , col] tut row* 1 , col+1 ] ; 
psiA  :=  (dely  *  dels  /  aZero)  *  (sum  /  4); 


This  copy  made  at  U.S.  Government  expense. 
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end: 


function  psiB:  real; 

yar 

topSum  real ; 

row,  col:  integer; 
begin 

topSum  : =  0 ; 
row  : =  1 ; 

for  col  : a  1  to  (nCols  -  1)  do 

topSum  :  =  topSum  ♦  u[row,col]  +  u[row, col+1] ; 
psiB  : =  (dely  /  bZero)  *  topSum  /  2 
tuad; 


function  heightVar (height:  real;  col:  integer):  real; 
var 

numerator, denominator  :  real; 
begin 

numerator  : =  height  *  height  *  bZero  * 

(1-  (bZeroPrime  *  aZero)  /  (bZero  ebZero  *  3)); 
numerator  :  =  numerator  *  u[l,col]; 
denominator  :s  2  *  aZero  *  (psiB  -  psiA) ; 
heightVar  : =  numerator  /  denominator 
end: 


procedure  surf aceVariat ion; 
begin 

for  col  : =  1  to  nCols  do 

begin  N 

height [col]  :s  (heightO  ♦ 

heightVar ( heightO , col ) ) /heightO ; 

writeln(col: 10,  height [col] : 10: 3)  > 

find  ; 

find;  r 

c 

procedure  compute; 
begin 

iterCount  : =  iterCount  +  1 ; 
topRow; 

for  row  :  =  2  to  (nRows-1)  dfi  midBow(row); 
bottomRow; 
normalise; 

writeln( iterCount: 10 ,  newSum:20:4,  newCorner : 10 : 3, 
convergeTest: 15: 6,  omega: 10: 3) 

And; 


function  finished:  boolean; 
begin 

convergeTest  : =  abs((newSum  -  oldSum)/newSum) ; 
converged  : =  convergeTest  <  tolerance; 
finished  :=  (converged  and  (iterCount  >=  minlterations ) ) 
or  (iterCount  >=  maxlterations) 

find; 


procedure  saveReaults; 
XAX 

widthStr:  str3; 
thetaStr:  str2; 
stepStr:  str3; 
heightStr :  str3; 


U.S.  (ioven 


r 
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s* 


filename:  stringri21 : 
savaFile:  text; 

procedure  f ill2(var  xStr:  atr2); 
begin 

while  length(xStr)  <  2  dfi  xStr  :=  ’0’  ♦  xStr 
sod; 

procedure  f 1113(var  xStr:  atr3); 
begin 

while  length(xStr)  <  3  dfi  xStr  :  =  ’O’  +  xStr 
and; 
begin 

str(width*10:0:0,  widthStr); 
atr (theta: 0:0,  thetaStr); 
at r ( sSteps : 0 ,  stepStr ) ; 
atr (heightO* 1000:0:0,  heightStr) ; 
f 1112 (thetaStr) ; 
fill3(wldthStr); 
f 1113 (stepStr) ; 
f 1113(helghtStr) ; 

fileName  :  =  thetaStr  ♦  widthStr  +  stepStr  ♦  ♦  heightStr; 

assign (aavaFl la,  fileName); 
rewrite (savaFi la) ; 

for  Index  :=  1  to  5  do  writaln (savaFi la ) ; 
wrltalnC savaFi la,  'Data  file:  ' ,  fileName); 

writeln( savaFile) ;  _  • 

writeln(saveFila,  'Depth  =  ',  dapth:0:l,  '  Bottom  width  =  ',  width: 0:1, 

’  Theta  -  ',  theta: 0:1,  '  Top  width  -  ',  topWidth'-0 : 1 ) ; 
writeln( saveFila ,  'delta  Z  =  dalZ:0:2,  '  delta  Y  r  ',  dalY:0:2, 

'  •  Row*  =  ',  nRows:0,  '  •  Colusms  =  ',  nCols:0); 

writaln (savaFile,  'Owaga  =  oasga:0:3, 

Convergence  tolerance  s',  tolerance); 
wrlteln( savaFile,  'Sub  over  matrix  =  new8ua:0:4, 

'  psiA  S  ',  psiA: 0:4,  '  psiB  =  ',  psiB:0:4); 

writaln (savaFile) ; 

write (savaFile,  iterCount:0,  '  iterations  completed;  computations  ’); 

If  converged  then  writeln(save?ile,  'converged.') 

else  writaln (savaFile,  'did  not  eonverge.'); 
writaln (savaFile,  'Relative  error  in  aatrix  sub  s  ',  converge Test ) ; 
writaln (savaFile,  duration:0:0,  '  seconds  for  iterations.'); 
writeln( savaFile) ; 
writaln (savaFile) ; 

writaln (savaFile,  'Transverse  surface  profile:'); 
write In ( savaFi le ) ; 

writeln( savaFile,  'Column': 12,  ’Y':7,  ’ Aaplitude' :  13) ; 
write In (savaFile) ; 

tSUL  col  '•  =  1  to  nCols  dfi  writaln  (saveFila,  col:  10, 

delY  *  (ool  -  1 ) : 10: 2,  height[col) : 10: 3) ; 

wr italn ( savaFile ) ; 
writeln( savaFile ) ; 

writeln( savaFile,  'Solution  matrix  (rotated  —  1,1  is  at  upper  right):’); 
writeln( savaFile) ; 
writaln ( savaFi le ) ; 
write (savaFile,  'ROW'); 

for  row  : =  nRows  downto  1  dfi  write( saveFila,  row: 8); 

writeln( savaFile) ; 

writaln (savaFile,  'COL'); 

writeln( saveFila) ; 

for  col  ::  1  is  nCols  do 

bSAiD 

write( savaFile,  col: 2,  '  '); 


